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This paper deals with extinction of an isolated population caused by intrinsic noise. We model the 
population dynamics in a "refuge" as a Markov process which involves births and deaths on discrete 
lattice sites and random migrations between neighboring sites. In extinction scenario I the zero 
population size is a repelling fixed point of the on-site deterministic dynamics. In extinction scenario 
II the zero population size is an attracting fixed point, corresponding to what is known in ecology 
as Allee effect. Assuming a large population size, we develop WKB (Wentzel-Kramers-Brillouin) 
approximation to the master equation. The resulting Hamilton's equations encode the most probable 
path of the population toward extinction and the mean time to extinction. In the fast-migration 
limit these equations coincide, up to a canonical transformation, with those obtained, in a different 
way, by Elgart and Kamenev (2004). We classify possible regimes of population extinction with 
and without an Allee effect and for different types of refuge and solve several examples analytically 
and numerically. For a very strong Allee effect the extinction problem can be mapped into the 
over-damped limit of theory of homogeneous nucleation due to Langer (1969). In this regime, and 
for very long systems, we predict an optimal refuge size that maximizes the mean time to extinction. 
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I. INTRODUCTION 

Every isolated population ultimately goes extinct. 
This happens, even in the absence of adverse environmen- 
tal variations, because of the discreteness of the individ- 
uals and random character of birth and death processes. 
Extinction risk is a major negative factor in viability of 
small populations 0, Hf, whereas extinction of diseases 
0, Q is usually beneficial. 

Extinction of a large population because of the intrin- 
sic noise demands an unusually large fluctuation: a rare 
sequence of random events when deaths dominate over 
births. Evaluating the role of rare large fluctuations in 
far-from-equilibrium systems is hard, and so population 
extinction, caused by intrinsic noise and environmental 
variations, has attracted much interest from physicists 
With a few exceptions 0, @, these studies as- 
sumed well-mixed populations, when spatial degrees of 
freedom are irrelevant. It has been known, however, 
since the classical paper of Skcllam j24[, that migration 
of individuals plays a crucial role in a host of natural 
environments of interest to population biology and epi- 
demiology [2o| . and in other applications. An important 
step forward in quantifying the extinction risk of spa- 
tially distributed populations was made by Elgart and 
Kamenev [(|. They considered a population on a dis- 
crete lattice that models a refuge of a large but finite 
size. The population undergoes on-site birth-death pro- 
cesses and migration of individuals between neighboring 
sites. Beyond the refuge the conditions are so harsh that 
they can be modeled by an infinite death rate. Elgart 
and Kamenev transformed the master equation for the 
evolution of a multi-variate probability distribution of 
the population size into an effective continuous classical 



mechanics by applying a time-dependent WKB (Wentzel- 
Kramers-Brillouin) approximation that uses the typical 
on-site population size K in the long-lived state of the 
population as a large parameter. The time-dependent 
WKB method yields a Hamiltonian functional and the 
corresponding Hamilton's equations - partial differential 
equations for an effective momentum p (coming from the 
probability generating function) and a conjugate coordi- 
nate q (that, in the deterministic limit, coincides with 
the populations size) . Both p, and q depend on the con- 
tinuous spatial coordinates x and time t. The extinction 
rate is determined by the classical action calculated along 
a special trajectory in the (infinite-dimensional) phase 
space q(x), p(x) of the system 

The present paper also deals with extinction of 
spatially-distributed populations caused by intrinsic 
noise. We suggest an approach that is closely related 
to that of Elgart and Kamenev Q , but also differs from 
it in a number of ways. First, in addition to scenario I 
of extinction, considered already in Rcf. Q, we also ad- 
dress scenario II. In scenario I the zero population size 
is a repelling fixed point of the on-site deterministic dy- 
namics. In scenario II it is an attracting fixed point, 
corresponding to what is known in ecology as Allee effect 
[26j. The results in these two extinction scenarios turn 
out to be quite different. Second, we derive the WKB 
equations systematically from the master equation for 
the multi-variate probability distribution. This deriva- 
tion shows that a continuous description in space is only 
valid when the migration rate between the neighboring 
sites greatly exceeds the on-site process rates. Third, by 
focusing on the long-lived quasi-stationary distribution 
of the population size, we formulate a stationary WKB 
theory in terms of the population size (treated as a "co- 
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ordinate") and its conjugate momentum. Fourth, an im- 
portant attribute of this WKB theory is spatial bound- 
ary conditions for WKB momentum p(x, t). We derive 
these boundary conditions, thus correcting an omission 
in Ref . Q . Fifth, using the WKB theory, we establish im- 
portant general properties of the most probable path of 
the population to extinction. We show that, in scenario 

I, the mean time to extinction (MTE) is determined by 
a hctcroclinic trajectory between two fixed points in the 
(infinite-dimensional) functional phase space of the sys- 
tem. The first fixed point corresponds to the long-lived 
quasi-stationary distribution of the population size. The 
second one corresponds to a zero-population-size state 
with a nontrivial momentum profile. In scenario II we 
only have results in the limit of a very strong Allee effect: 
close to a characteristic bifurcation of the system. Here 
again we obtain the solution of the problem in terms of a 
heteroclinic connection: between the fixed point, corre- 
sponding to the long-lived quasi-stationary distribution, 
and a fixed point describing the "critical nucleus" . In this 
limit the population extinction problem turns out to be 
completely integrable, similarly to the integrability of the 
problem of population explosion close to the saddle-node 
bifurcation We explain this integrability by estab- 
lishing a direct connection between this problem and the 
over-damped limit of theory of homogeneous nuclcation 
due to Langer [27[ ■ We consider different types of refuge, 
determined by the conditions at the refuge boundaries 
and illustrate our results by solving, analytically and nu- 
merically, three particular population models. In most 
of this paper we deal with refuges whose spatial sizes are 
not exponentially large in parameter K 3> 1. An ex- 
ception is section |V] B where extinction of populations 
residing in very large refuges is considered (again, for a 
very strong Allee effect). Surprisingly, we find here an 
exponentially large reduction in the MTE and predict an 
optimal refuge size that maximizes the MTE. 

The remainder of the paper is organized as follows. 
Section [H] includes important preliminaries which are 
used in the subsequent sections. It gives an overview of 
deterministic theory of population dynamics in a refuge: 
with and without Alice effect, and for different spatial 
boundary conditions. It also discusses, on a qualitative 
level, how the noise-driven population extinction is ex- 
pected to occur in different cases. Section ITU presents a 
stochastic theory of the population dynamics in a refuge. 
Here we introduce the master equation, focus on the 
quasi-stationary multi-variate distribution of the popu- 
lation sizes and on the MTE, and formulate a WKB the- 
ory aimed at evaluating these quantities. Sections IIVI 
and [V] analyze population extinction in scenarios I and 

II, respectively. Here we consider two specific birth-death 
models in the region of parameters close to their charac- 
teristic bifurcations. In this way we achieve some gener- 
ality, as the reduced equations, in each of the two cases, 
describe a broad class of population models. We also 
revisit, in scction llVBl an additional model problem, ex- 
hibiting extinction scenario I. Extinction of populations 



residing in exponentially large refuges is considered, for 
a very strong Allee effect, in section [V] The results are 
discussed, along with some possible generalizations and 
unresolved problems, in section [VTl 

II. DETERMINISTIC EQUATIONS AND 
POPULATION EXTINCTION SCENARIOS 

A. General 

Consider a single population residing in a refuge by 
which we mean a one-dimensional lattice of N ^> 1 sites 
(or habitat patches) labeled by index i = 1,2, ...,N. 
The population size rii at each site varies in time as a 
result of two types of Markov processes. The first set 
of processes involves a local, on-site stochastic dynam- 
ics of birth-death type, with birth and death rates A(rij) 
and n(rii), respectively, where fi(Q) = 0. As there is 
no creation of new individuals "from vacuum" , one has 
A(0) = 0. The second process is random and indepen- 
dent migration of each individual between neighboring 
sites with migration rate coefficient Dq- What happens 
at the edges of the refuge, i = 1 and i = N, needs to 
be specified separately; we will deal with this issue a bit 
later. 

Assuming m S> 1, one can attempt to neglect fluctua- 
tions and describe the population dynamics by determin- 
istic rate equations 

hi = A(rij) - ju(ni) + Do(n,-i + n l+1 - 2rij) . (1) 

Established populations are described, in the determin- 
istic limit, by stable steady-state solutions n, of this set 
of TV coupled equations. According to Eq. (fTJ), an estab- 
lished population would persist forever. The stochastic 
picture is markedly different. An unusual sequence of 
births and (predominantly) deaths can bring the popula- 
tion to the absorbing state (rii = 0, ri2 = 0, . . . , njv = 0) 
corresponding to extinction occurring everywhere. This 
ultimately happens with probability one. 

Before dealing with the stochastic problem, however, 
let us dwell some more on deterministic rate equations (TTJ) 
and their predictions. Let the characteristic population 
size on a single site, predicted by a steady-state deter- 
ministic solution scales as K ^> 1 . This implies 0, H3| 
that, in the leading order of K, one can represent the 
birth and death rates as 

A(rij) = jU A'A(ft) and fi(rn) = fj l0 Kp(q i ), (2) 

where qi = rii/K is the rescaled population size at site i, 
A (ft) ~ P-{qi) ~ 1) and fiQ is a characteristic rate coeffi- 
cient. Now Eq. (TTJ can be rewritten as 

ft = Mo/(ft) + A)(ft-i + ft+i - 2ft) , (3) 

where /(ft) = A(ft) — £t(ft) is the rescaled birth-death 
rate function. 
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With no migration, Do = 0, the on-site deterministic 
dynamics is determined by the equation q = Haf(q). One 
fixed point of this equation is q = 0, and there are two 
major cases determined by the sign of derivative /'(<?) 
at g = 0. For /'(0) > (scenario I) the fixed point 
q = is repelling, and the on-site population size, in the 
absence of migration, flows to an attracting fixed point 
q = qi > that describes an established population. 
One example of scenario I is the well known SIS model 
of epidemiology [28[ for which X(n) = X n(K — n) and 
/x(n) = fiQn. Here X(q) = Roq(l — q), p.(q) = q, and 
/(g) = q(Ro — 1 — Roq), where Rq = XqK/^lq is the basic 
reproduction number. At Rq > 1 q = is a repelling 
point of equation q = /io /(<?), whereas q = qi = 1 — 1/Ro 
is an attracting point. 

In scenario II one has /'(0) < 0. Here fixed point q = 
is attracting, and the population gets established, at an- 
other attracting fixed point q = q2 , only if the initial pop- 
ulation size exceeds a threshold: a repelling fixed point q\ 
such that < <7i < q-2- Scenario II accounts, in a simpli- 
fied way, for a host of Allee effects [26|. As an example of 
scenario II we will consider the following three reactions: 
A -> 0, 2A -> 3A and 3A -> 2A with rate coefficients 
/io, Ao and a , respectively 0, [H|. Here X(q) = 2q 2 /^ 
and p,(q) = q(l + q 2 /^), where K = 3A /(2cro) and 
7 = 8/!oCo/(3Ao). At S 2 = 1 — 7 > the system exhibits 
bistability. Here the zeros of the rescaled birth-death rate 
function 

f{q) = --q{q-qi){q-q2) (4) 
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describe two attracting fixed points, and q2 = 1 + 6, and 
a repelling fixed point qi = 1 — 6 such that < q\ < q^. 

Now let us reintroduce deterministic migration and as- 
sume that it is much faster than the on-site population 
dynamics: D a 3> (J-q (the criterion can become less re- 
strictive close to characteristic bifurcations of the on-site 
population models, see sections [IV Al and |V|) . In this case 
one can use a continuous spatial coordinate x instead 
of the discrete index i and replace the discrete Lapla- 
cian in Eq. ([3]) by the continuous one. This brings about 
reaction-diffusion equation 

d t q = ^f(q) + Dd 2 x q, (5) 

where D = D$h 2 is the diffusion constant, and h is the 
lattice spacing. The system size becomes L = Nh. Equa- 
tion ([S]) ; which has been the subject of numerous stud- 
ies [25|, |3(|, should be supplemented by spatial bound- 
ary conditions. We will separately consider periodic, 
q(x + L) = q(x), and zero, q(0) = q(L) = 0, bound- 
ary conditions. In the discrete version of the problem, 
the zero boundary conditions correspond, up to small 
corrections (see Appendix A), to absorbing boundaries 
at sites i = 1 and i = N. The absorbing boundaries 
model, for example, extremely harsh conditions outside 
of the refuge [f| [24 |. Results for still another type of 
boundaries - reflecting walls at x = and x = L - can 



be easily obtained from the results for periodic boundary 
conditions. 

Spatial profiles of established populations are de- 
scribed, in the deterministic theory, by stable steady- 
state solutions q = q(x) > of Eq. ([5]). They satisfy 
ordinary differential equation 

Dq"(x)+fx f(q) = (6) 

subject to the chosen spatial boundary conditions. The 
first integral of this equation, 

^-{q'f + V{q) = const, (7) 

with effective potential V(q) ~ /(£) makes the 
problem soluble in quadratures and yields a phase por- 
trait of the steady states on the plane (q,q')- No- 
tably, reaction-diffusion Eq. ([5]) is a gradient flow, dtq = 
—5J-/5q, where 

F[q{x, t)} = f dx [-fi V(q) + (1/2) D{d x q) 2 } . (8) 
Jo 

Therefore, it describes a deterministic flow towards a 
minimum of the Ginzburg-Landau free energy F[q]. This 
property helps identify linearly stable and unstable in- 
dependent solutions, as they correspond to local min- 
ima and maxima of T[q] , respectively [13 • Furthermore, 
it yields a simple selection rule in cases when, at fixed 
L, there are multiple solutions of Eq. © with periodic 
boundary conditions: the solution with the maximum 
spatial period (equal to L) is selected when starting from 
a generic initial condition [3~i| . 




q 



FIG. 1: (color online) Effective potential V(q) and phase por- 
trait (q, q') for steady-state solutions of Eq. © in scenario I 
(no Allee effect). 



B. Scenario I 

What is the steady state in scenario I, as exemplified 
by the spatio-temporal SIS model? Figure Q] shows ef- 
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FIG. 2: (color online) Steady-state solutions of Eq. ((5]) in 
scenario I (no Allee effect) for periodic (a) and zero (b and c) 
boundary conditions in space. The arrows indicate extinction 
transitions driven by rare large fluctuations. L = 1.1L C (b) 
and 4L C (c). R = 2, so L c = 7r(D//i ) 1/2 , see Eq. ©. 



fective potential V(q) = (Rq - l)q 2 /2 - Roq 3 /3 and the 
resulting phase portrait (q, q') at Rq > 1. The only non- 
trivial steady-state solution, obeying periodic boundary 
conditions, is the x-indepcndcnt solution q = q±, depicted 
in Fig. Introducing intrinsic noise, we will see that 
the most probable path of this population to extinction 
is such that the population size drops to zero uniformly 
on the whole interval < x < L. For large systems, the 
MTE is very long in this case, being exponentially large 
in KL/h = KN. 

For the zero boundary conditions, an x-dependent 
steady state corresponds to a phase trajectory inside the 
scparatrix in Fig. [TJ Such steady states, depicted in 
Figs. [2] b and c, exist only if the system size L is larger 
than critical size 
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This quantity can be obtained from Eq. ^ linearized 
around q = 0. At L < L c there is only trivial solu- 
tion: no established population. The ^-dependent solu- 
tion emerges, at L = L Cl via a transcritical bifurcation. 
At L 3> L c the population size is close to q\ everywhere 
except in boundary layers, with thickness of order of L c , 
at x = and x — L. At L > L c the most probable path to 
noise-driven extinction for the zero boundary conditions 
is such that the population size drops to zero uniformly 



on the whole interval < x < L. As we will see in sec- 
tion |IV] (see also Ref. [|), for L > L c the MTE is again 
exponentially long in parameter KN. It becomes much 
shorter as L approaches L c , see section HVl 



C. Scenario II: Allee effect 

Now consider scenario II, on the example of three re- 
actions A -> and 2 A <=^ 3 A. At < 7 < 1, that is 
< S < 1, the effective potential, 
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has two maxima: at q = and q = q2, and which of the 
steady-state solutions q = and q = o 2 "wins" depends 
on which of the maxima is higher [25l . l30l | . 



1. Strong Allee effect 

Figures |U and H] illustrate tne case of V(q 2 ) < V(0): 
a strong Allee effect. In our example this occurs at 
8/9 < 7 < 1, or < S < 1/3. For periodic bound- 
ary conditions, the only linearly stable nontrivial steady- 
state solution is the x-independent solution q = qi- A 
sufficiently large perturbation, however, triggers a deter- 
ministic transition from q = q2 to the trivial solution 
q = that is also linearly stable. An important attribute 
of this metastability is presence of the "critical nucleus" : 
an a;-dependent solution q c (x) of Eq. ^ that is linearly 
unstable under the dynamics of Eq. ([5]). A small per- 
turbation around the critical nucleus brings the system 
either to q — or to q = q 2 . The critical nucleus is se- 
lected by the system size L and corresponds to a phase 
trajectory inside the internal separatrix shown in Fig. [3}}. 
The critical nucleus exists only for L > L Cl where 
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as can be obtained from Eq. ((6j), linearized around 
q = qi, with periodic boundary conditions. At L 3> L c 
the critical nucleus coincides with the internal separatrix 
in Fig. For f(q) from Eq. ^ (a cubic polynomial), 
the critical nucleus can be found analytically, in terms 
of elliptic functions, by integrating the first-order equa- 
tion and choosing the solution q(x) with period equal 
to the system size L. A more practical way is to solve 
Eq. (|6]) numerically, by shooting. One solves numerically 
an initial-value problem for Eq. ^ starting, at x = 0, 
from some q(0) £ (gi,^) and q'(0) = 0. Parameter q(0) 
is varied until the numerical solution exhibits a single 
full-period oscillation, so that q(L) ~ q(0) and q'(L) ~ 0. 
Figure U shows the critical nuclei for two different values 
of L > L c . Note that a critical nucleus corresponds to a 
local maximum of free energy (|73p (30j . 
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The presence of a critical nucleus in the deterministic 
theory plays a pivotal role in the noise-driven extinction 
of an established population exhibiting a strong Alice 
effect. Indeed, a large fluctuation of the size of stochas- 
tic population residing around q = q 2 can create critical 
nucleus q c (x). The further population dynamics toward 
extinction proceeds "downhill" , that is essentially deter- 
ministically. What happens at L ^> L c , sec Fig. @J>, is 
intuitively clear, and will be supported by our quantita- 
tive results in section [V] Here the rate of noise- induced 
creation of the critical nucleus is exponentially small in 
K but independent of L (unless L is exponentially large 
in K). Once having passed the critical nucleus, the so- 
lution q(x, t) of Eq. dU develops, on a time scale ~ ^ 1 , 
two outgoing deterministic "extinction fronts" . In our 
example of three reactions the deterministic front solu- 
tion can be found analytically [25l l30l ] . The extinction 
fronts propagate with speed 
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and drive the whole population to extinction on a time 
scale ~ Lf(fioD) 1 / 2 . Therefore, unless the system size L 
is exponentially large in K, it is the creation of a single 
critical nucleus that serves as the extinction bottleneck. 
That is, the MTE is determined here by the mean cre- 
ation time of the critical nucleus. This quantity does not 
include an exponential dependence on the system size 
L and is therefore much shorter than in scenario I. Now, 
what happens when L is above L c but close to it? We will 
show that here too the most probable path to extinction 
corresponds to a large fluctuation bringing the popula- 
tion from q = q 2 to critical nucleus q c {x), see Fig.@Jt, and 
not to the x- independent unstable state q = qi . 
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FIG. 4: (color online) Linearly stable states q — qi and q = 
0, linearly unstable state q — qi and critical nucleus q c (x) 
for a strong Allee effect and periodic boundary conditions. 
The system size L — 1.04L C (a) and 3.7 L c (b), where L c is 
defined in Eq. (|11|) . The arrows indicate transitions, driven 
by rare large fluctuations and leading to a rapid extinction. 
Parameter 7 = 24/25, so 6 = 1/5, and L c = 2n(3D/ A*o) 1/2 • 



For a strong Allee effect and zero boundary conditions, 
there is only (linearly stable) trivial steady state q = 0: 
no established population. 




FIG. 3: (color online) Effective potential V(q) and phase por- 
trait (?,<?') for steady-state solutions of Eq. ((5} for a strong 
Allee effect, V(0) > V(q 2 ). 



FIG. 5: (color online) Effective potential V(q) and phase por- 
trait (q, q') for steady-state solutions of Eq. ((5| for a weak 
Allee effect, V(0) < V{q 2 ). 
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FIG. 6: (color online) Linearly stable steady states q = q2 and 
q = 0, linearly unstable state q = qi and critical nucleus q = 
q c (x) for a weak Allee effect and periodic boundary conditions. 
The system size L = 1.4L C (a) and 3.25L C (b), where L c is 
defined in Eq. (|f II) . The arrows indicate extinction transitions 
driven by rare large fluctuations. 7 = 3/4, so 5 = 1/2 and 
L c = 7r(6D/ M o) 1/2 . 
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FIG. 7: (color online) The s-dependent linearly stable steady 
state (the upper curve), trivial stable state q = 0, linearly 
unstable state q — qi and critical nucleus q = q c (x) for a weak 
Allee effect and zero boundary conditions for L = 1.1L C (a) 
and 2.52 (b). The arrows indicate the extinction transitions 
driven by rare large fluctuations. 7 = 3/4, so S = 1/2, and 
L c ~ 6.026(D/ Mo ) 1/2 . 



2. Weak Allee effect 

For a weak Alice effect one has V(q2) > V(0), as illus- 
trated in Figs. [5] and [5] In our example of three reactions 
this case corresponds to < 7 < 8/9, or 1/3 < 5 < 1. 
For periodic boundary conditions there are two linearly 
stable ^-independent steady states, q = q 2 and q = 0, 
and the linearly unstable x- independent state q = q\. 
There is also critical nucleus q = q c {x), described by a 
phase trajectory located inside the internal separatrix in 
Fig. [5)3; it is selected by the system size L. The criti- 
cal nucleus exists when L > L c , where L c is given by 
Eq. (|TT|) . At L ^> L c the critical nucleus is described 
by the internal separatrix of Fig. [SJ). Here the popula- 
tion size, corresponding to the critical nucleus q c (x), is 
close to zero everywhere except in a narrow region with 
thickness ~ L c . What is the most probable path of the 
population toward noise-driven extinction? Here one has 
to choose between two paths. In the first path the pop- 
ulation size goes down from q = q2 to the x-independent 
unstable state q = qi on the whole interval < x < L 
and then continues to fall, almost deterministically, to 
zero. In the second path the population size goes down 
from q = qi to the critical nucleus and then, almost de- 
terministically, to zero. For L 3> L c the MTE involves, 
for each of the two options, an exponential dependence 
on L, so it can be very long. 

For the zero boundary conditions there are two lin- 



early stable steady states: an x-dependent state and the 
trivial state q = 0. There is also critical nucleus: an 
x-dependent unstable steady state. These solutions are 
depicted in Fig. [7J a and b. Each of the x-dependent solu- 
tions is selected by the system size L and described by a 
phase trajectory located between the two separatrices in 
Fig. [SJd. Among them there is a limiting phase trajectory 
such that the stable steady state, at given L, corresponds 
to a phase trajectory located between the limiting phase 
trajectory and the external separatrix. In its turn, the 
critical nucleus, for the same L, corresponds to a phase 
trajectory that lies between the limiting phase trajec- 
tory and the internal separatrix. The x-dependent solu- 
tions, both stable and unstable, exist when the system 
size L is larger than a critical size L c [which is different 
from L c given by Eq. (fTTj) ]. The critical size L c scales 
as (D/^o) 1 ^ 2 and also depends on S. For L >• L c the 
linearly stable steady state corresponds to the external 
separatrix and is therefore close to qi everywhere except 
in boundary layers with thickness ~ L c at x = and 
x = L. In its turn, the critical nucleus corresponds, at 
L 3> L c , to the internal separatrix and therefore coincides 
with the critical nucleus obtained for periodic boundary 
conditions. At L = L c the stable and unstable solutions 
merge. At L < L c there is only trivial steady state q = 
which is linearly stable. The most probable path to ex- 
tinction at L > L c corresponds to a large fluctuation that 
brings the population size from the stable state down to 
the critical nucleus, see Fig. [7] a and b. 
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III. MASTER EQUATION AND WKB 
APPROXIMATION 



A. Governing equations 

Now let us return to the discrete-lattice model and de- 
scribe stochastic dynamics of the population. This can be 
done in terms of evolution of the multivariate probabil- 



J 



ity distribution P(n,t) — P(ni, ri2, . ■ . , t) = P(n, n,, t), 
where i = 1,2, . . . , N , and n denotes the vector of all n's 
not explicitly written, see e. g. Ref. [H]. This probability 
distribution is assumed to be identically zero if any of 
is negative. For the continuous-time Markov processes 
of birth, death and migration, the master equation for 
P(n, t) has the following form: 



i=i 



d t P{n, t) = J2 { A ("« - 1)-P(n, nt -l,t) + ti{m + l)P(n, n, +l,t)- [\(m) + p(ni)]P(n, *)} 
' i 

N 

,^{(ni_i + l)P(n,ni_i + l,m -!,*) + (n*+i + 1)P(^,th - l,n i+ i + l,t)- 2riiP(n,t)} . (13) 



This equation holds as it is for a periodic lattice with 
period N. For absorbing boundaries the migration terms 
i = l and i = N arc slightly different, see Appendix A. Of 
a primary interest for us is the instantaneous extinction 
rate, or extinction probability flux: 



d t p(o,t) = n(i)[p(i,o,...,o,t) 
+ P(0,0,...,M)]- 



P(0,1, 



,0,t) + ... 
(14) 



We will continue to assume that K 3> 1. Furthermore, 
we will assume in most of the paper (except in section 
|V]B) that the system size is not too large: not exponen- 
tially large in K. In this case, extinction of an established 
population proceeds, in the probabilistic language, as fol- 
lows. During the relatively short relaxation time t r , de- 
termined by the deterministic rate equation (TTJ), the sys- 
tem approaches a quasi-stationary state, where P(n, t) 
is sharply peaked at the relevant steady-state solution of 
Eq. (JTJ) . At t ^> t r the quasi-stationary probability slowly 
decays in time, 



P(n,i) ~ 7r(n)e 



-*/T e 



0,1,2, 



(15) 



r 



except for n = that corresponds to a complete extinc- 
tion. The decay rate 1/T e is the lowest positive eigen- 
value of the time-dependent master equation (|13[) . This 
eigenvalue is special: it turns out to be exponentially 
small with respect to A' > 1 j33|. The probability 
of complete extinction P(0, t) = P(0, 0, . . . , 0, t) slowly 
grows in time: 



P(0,t) ~ 1-e 



-t/T e 



(16) 



In this regime the decay time T e is equal to the MTE, 
whereas the probability distribution of extinction times 
is an exponential distribution with mean T e , see e.g. Ref. 
fl~0| . Using Eqs. (JTSJ) and (fill), wc can rewrite Eq. (fT3|) 
as a linear eigenvalue problem for the quasi-stationary 
distribution 7r(n): 



^ |A(?ii - 1) 7r(n, rii - 1) + n(m + l)7r(n, n.; + 1) - [\(rn) + //(n,)] 7r(n)| 

i=i 

N 

+D ^2{(ni-i + 1) n(n,n,i-i +1,^-1) + (n i+1 + l)n(h,n t - l,n i+1 + 1) - 2n i7 r(n)| = -Att(ii) , (17) 



(except for n = 0) for the lowest positive eigenvalue A = relation 
1/T e . Once 7r(n) is determined, A can be found from 

A = ^(l)[7r(l,0,...,0) + 7r(0,l,...,0) + ... 

+ tt(0,0,...,1)] (18) 
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following from Eqs. (|14)) - (fT6|) . 

For K 3> 1 and n.; > 1 we can treat qi = rii/K as 
continuous quantities and solve Eq. (fT7| in WKB ap- 
proximation which generalizes to s patial pop ulations the 
stationary WKB method [|,[il,[i3M[lI[2J,[3l|3i| pre- 
viously employed for well-mixed populations. The WKB 
ansatz is 



7r(n) = cxp [-KS(q)} . 



(19) 



Our goals are to accurately evaluate the leading-order 
contribution to ln(/x r e ) and to find the most probable 
path of the population to extinction. We plug Eqs. © 
and (fT9|) in Eq. (fT3|) and neglect term — A7r(n) which is 
expected to be exponentially small in K 1. In the 
leading order in 1/K this procedure yields a stationary 
Hamilton- Jacobi equation 



H(q,d^S) =0 



(20) 



with an effective classical Hamiltonian with N degrees of 
freedom, 



N 



#(q,P) =^o£ [Hu) (e Pl - 1) + P(qi) (e- pi - l)] 

JV 

+D J2h-i (e*-*- 1 - l)+%+i (e p '- p ' +1 - l)] ,(21) 



where p.; = d qi S. This lattice Hamiltonian, and corre- 
sponding Hamilton's equations - a set of 2N ordinary 
differential equations for qi(t) and Piit) - is a proper 
framework for dealing with population extinction for any 
relation between the migration rate coefficient Dq and the 
characteristic rate coefficient po of the on-site dynamics 

In the following we will only consider the limit when, 
as in section UH migration between the neighboring sites 
is much faster than the on-site population dynamics: 
Do 3> Mo (the criterion becomes softer close to bifurca- 
tions of the on-site models, see sections llV Al and W} . In 
this regime the quasi-stationary distribution 7r(n) and, as 
a consequence, the classical action S'(q) are slowly vary- 
ing functions of n and q, respectively. This implies that 
the difference between the momenta pi on neighboring 
sites is much smaller than unity. Taylor-expanding the 
migration term H m of Hamiltonian (|2ip (the term pro- 
portional to Do) up to second order, we obtain 

N 

H m (q, p) = D [- (Qi - (Pi - Pi-i) 

i=l 

+ \{qi + qi-i)(jPi-Pi-i?]- (22) 

The slow variation of qi and Pi with i calls for a con- 
tinuous description. We introduce a continuous spatial 
coordinate x instead of the discrete index i and arrive at 



an effective continuum classical mechanics. The Hamil- 
tonian functional is 

H[q(x,t),p(x,t)] = ~ [ dxw, (23) 
" Jo 



with density 



w = H a {q,p) - D d x qd x p - q (d x p) (24) 

and on-site Hamiltonian 

H (q,p) = Mo [%) (e p - 1) + TM) (e- p - l)] . (25) 

Note the presence of two diffusion terms inside the square 
brackets in Eq. (f24|) . The first term describes determin- 
istic diffusion, the second one describes fluctuations of 
diffusion. Hamiltonian, related to Eq. (|2"3")l by canon- 
ical transformation Q = qe~ p , V = e p , was obtained 
by Elgart and Kamenev [6j who employed the proba- 
bility generating function in conjunction with a time- 
dependent WKB theory. Note that the two diffusion 
terms in Eq. (|24|) add up to —D d x Qd x V in canonical 
variables Q and V . This simplification, and the some- 
what simpler form of the on-site Hamiltonian, can be 
advantageous, see section ITVBl 

The Hamilton's equations of motion, 



O t q = fl- 
op 



d t p 



5H 



Ho [Mq)e p - fi(q)e- p ] 

D [&lq - 2d x (qd xP )] , (26) 



d d 2 xP + (d xP y 



(27) 



are partial differential equations for continuous variables 
q(x,t) and p{x,t) = h5S/Sq (37} • Note that, for the 
purpose of solving stationary Hamilton- Jacobi equation 
(|20| . time as appears in Hamilton's Eqs. (|26| and (f27f is 
merely a way of parametrizing phase space trajectories. 
It is not necessarily related to the original time entering 
Eqs. (|13p -([T6 |) for the evolution of probabilities. To re- 
mind the reader, Eqs. ([2T|) and (|22|) . as well as Eqs. (fl"3|) 
and (|17|) . are only valid for periodic systems; absorbing 
boundaries are considered in Appendix A. It is impor- 
tant, however, that continuous equations (f2"3"]) - (|2"7) are 
valid in the case of absorbing boundaries as well. 

For all types of spatial boundaries, continuous 
Eqs. (|26| and ([27]) must be complemented with spa- 
tial boundary conditions. This circumstance was left 
unattended in Ref. @. For periodic systems the spa- 
tial boundary conditions are or course q(0,t) = q(L,t) 
and p(0,t) = p(L,t). For reflecting boundaries they 
are also straightforward: d x q(0,t) = d x q(L,t) = and 
d x p(0,t) = d x p(L,t) = 0. The case of absorbing bound- 
aries is a bit more involved, and we derive the correspond- 
ing boundary conditions in Appendix A. Up to small cor- 
rections O(no/ Do) 1 / 2 <C 1, they turn out to be zero con- 
ditions both for the coordinate, and for the momentum: 
q(0, t) = q(L, t) = 0, and p(0, t) = p(L, t) = 0. 
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B. Activation trajectories 

Now let us return to Eqs. (fT9| and (|20j) that describe, 
in WKB approximation, the quasi-stationary distribu- 
tion 7r(n). This distribution is smooth and has its (Gaus- 
sian) maximum at q(x) — q s (x). Therefore, in order 
to find 7r(n), one needs to find a particular solution of 
Hamilton- Jacobi Eq. (f2"U)) such that its variational deriva- 
tive vanishes at q = q a (x): SS/Sq\ q ^^ = 0. Setting 
5 {q s (x)} = 0, we define 5 {q(x)} uniquely as the solu- 
tion of Eq. (|20p. Once this solution is known, one can 
use Eq. (fT8|) to evaluate the MTE up to prc-cxponential 
factors: 

ln( Mo T e ) ~ KS(0) . (28) 

In order to calculate 5(0) we will use Hamilton's equa- 
tions (f2l)]) and (j2~7| that describe trajectories in the func- 
tional phase space {q(x),p(x)}. As Hamiltonian (j2"5)l 
does not depend explicitly on time, it is a constant 
of motion. Furthermore, as Hamilton- Jacobi equation 
([2"0")l is stationary, we should only consider trajectories, 
for which this constant of motion - the total energy of 
the effective mechanical system - is zero. The simplest 
among zero-energy trajectories are deterministic, or re- 
laxation, trajectories: solutions of Eqs. (|2"o]l and (|2"7|) with 
p{x, t) = 0. Here Eq. (|2"o]l reduces to the deterministic 
reaction-diffusion equation (O, whereas Eq. (f27|) is satis- 
fied trivially. 

The quasi-stationary distribution 7r(n) is peaked at 
what we call fixed point A: (functional) fixed point 
q(x) = q 8 (x), p{x) = of Eqs. J26j and There- 
fore, the phase trajectory we are interested in for the 
purpose of calculating 5(0) should start, at t = — oo, 
at fixed point A. In both extinction scenarios I and II 
there are a stable manifold p{x) = 0, and an unstable 
manifold p(x) ^ 0, emanating from fixed point A, see 
Appendix B. In the discrete lattice formulation, each of 
these two manifolds is A^-dimcnsional and is embedded 
into zero-energy hyper-surface H{q(x),p(x)} = 0. 

For any phase trajectory that originates from fixed 
point A at t = — oo, we can write 

S{q(x,T)} = - [ dt [ p(x,t) d t q(x,t) dx . 

h J-oo Jo 

In view of Eq. (|28|i , we only need to consider phase tra- 
jectories that reach extinction hyper-plane q(x) = [so 
that q(x,t) vanishes at all x]. Relaxation trajectories, 
p = 0, that exit fixed point A, cannot reach the extinc- 
tion hyper-planc, so we need an activation trajectory, 
p / 0, for this purpose. For extinction scenario I, a 
crucial property of the activation trajectory can be es- 
tablished under quite general assumptions. The activa- 
tion trajectory must approach, at t — +oo, another fixed 
point which we call fixed point B. It involves q(x) = 0, 
see Fig. [21 and p(x) = p e {x): the non-trivial steady-state 
solution of Eq. (|2"T)) with q(x) = and proper spatial 



boundary conditions. That is, in scenario I the activa- 
tion trajectory must be a heteroclinic connection AB [or 
instanton, see Ref. [H] for a review on instantons] in 
functional phase space {q(x),p(x)}. The proof of this 
statement is presented in Appendix B2; it relies on the 
structure of the phase space of Eqs. (|26|) and (|27|) and, in 
particular, on the presence and linear stability properties 
of (zero-energy) fixed points of Eqs. (|2"rJ)) and (|2"T)l . 

For extinction scenario II the structure of the phase 
space is more complicated, and we cannot make an 
equally general statement about the properties of the 
activation trajectory, except that this trajectory must 
exit, at t = — oo, fixed point A and ultimately arrive 
at extinction hypcrplanc q = 0. We know much more, 
however, in the case of a very strong Alice effect, when 
the basin of attraction of the state q — qi in the de- 
terministic theory is small. Here the noise only needs 
to create the critical nucleus, see Section II CI. In the 
language of Eqs. (f2T)]) and (|27[) . the activation trajectory 
must approach, at t = +oo, fixed point D that involves 
q = q c {x) (the critical nucleus) and p(x) = 0. One can 
argue that, from there on the population flows to fixed 
point C (where q = p = 0) along a relaxation trajectory. 
The relaxation trajectory does not cost any action (un- 
less the system size L is exponentially large in K, see 
section |V]B). In this case the MTE can be identified, up 
to a pre-exponent, with the mean time of creation of the 
critical nucleus. Note that, in this limit, the activation 
trajectory is again a heteroclinic connection (AD). One 
can expect that, for a moderately strong Allee effect, the 
activation trajectory will still involve a critical nucleus 
and, therefore, represent a heteroclinic connection AD. 

Once the activation trajectory is found, we can obtain 
the MTE in the leading order of the WKB theory by 
calculating 5(0), entering Eq. (|2"5]l , along this trajectory. 
In scenario I 5(0) is the action along the heteroclinic 
connection AB. In the strong- Alice-effect limit of scenario 
II one has ln(/xoT e ) ~ KSq, where 5o is the action along 
the heteroclinic connection AD. In both cases we can 
write ln(^o^e) — KS, where 

1 f^ 1 f°° 
S=j dx dtp(x,t)d t q(x,t) . (29) 

If there are more than one heteroclinic connections be- 
tween the same pair of fixed points, and obeying the same 
boundary conditions in space, one should choose the con- 
nection which yields the minimum action. Similarl y t o 
non-spatial but multi-population systems [rl H2L [l3l I22I] . 
the minimum-action trajectory is the most probable path 
of the population on the way to extinction. Sections IIVI 
and [V] present three particular examples of determining 
the activation trajectories and evaluating the MTE. 

As we already mentioned, Hamilton's equations (|2"o]l 
and (|27p coincide, upon canonical transformation Q = 
qe~ p , V = e p , with those derived by Elgart and Kamenev 
in the framework of a time-dependent WKB approxima- 
tion Q. There are some differences, however, between 
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our and their formulations of the problem. These differ- 
ence involve boundary conditions: both in time, and in 
space. The differences in the boundary conditions in time 
appear already in the most basic, spatially-independent 
setting, so let us consider this setting first. 

The time-dependent WKB formulation of Ref. [f| pre- 
scribes, at t = 0, the initial population size, say q = q s , 
with an a priori unknown momentum. It also prescribes, 
at a (sufficiently large) final time t = T, momentum 
V = (or, in our variables, p = -co), with an a priori 
unknown population size. One needs to find the initial p 
and the final q from the condition that the action along 
the resulting trajectory is minimum. T is ultimately sent 
to infinity [6[. 

Our WKB formulation differs, first of all, in its pre- 
scription of the final state of the system. In view of 
Eq. (|2"5|) , we demand q — there. Furthermore, we know 
that the activation trajectory must be a heteroclinic con- 
nection AB (in scenario I) or AD (in scenario II). This 
involves a full knowledge of both q and p at the initial 
(t = — oo) and final (t = oo) points. Importantly, the 
final value of the momentum in this formulation is differ- 
ent from p = — oo, or V = demanded in Ref. jf|. 




FIG. 8: Zero-energy trajectories in spatially- inde pen dent set- 
tings for extinction scenario I (a) and II (b) [20l. |29|. Shown 
are fixed points A,B and C (a) and A,B,C and D (b). The 
arrows show stable and unstable manifolds of the correspond- 
ing fixed points. The activation trajectories AB (a) and AD 
(b) are accentuated by thicker lines. 

In spite of these differences, the two formulations yield, 
in the spatially-independent case, the same result for the 
MTE. This happens because of two reasons. First, the 
activation trajectory that emerges, at T — > oo, in the 
time-dependent formulation has zero energy, as in 
our quasi-stationary formulation. Second, the activation 
trajectory in the time-dependent formulation is, in gen- 
eral, not a heteroclinic connection. Rather, it consists 



of two (in scenario I) or even three (in scenario II) sep- 
arate heteroclinic connections [13, [H], see Fig. [H One 
of them coincides with trajectory AB or AD (for scenar- 
ios I or II, respectively) that the quasi-stationary theory 
predicts. The other segments go along either q = or 
p = lines and therefore do not contribute to the action. 
[One can even argue that the last segment ultimately 
reaches p = — oo, see Fig. [5J] An advantage of the quasi- 
stationary theory, especially in numerical calculations, is 
that the non-contributing segments of the trajectory are 
excluded from the start. 

For spatially-dependent systems the differences be- 
tween the two formulations may become irreconcilable. 
Consider, for example, scenario I in the case of absorbing 
boundaries. Here the x-independent momentum p(x) = 
— oo, postulated as the final state in Ref. @, does not 
obey the zero boundary conditions p(0,t) — p(L,t) = 0, 
and so it cannot possibly be a correct final state. 

IV. POPULATION EXTINCTION: SCENARIO I 

A. Universal limit 

Here we consider, as an example, the spatio-temporal 
SIS model. To render our results a broader relevance, 
we assume from the start that the basic reproduction 
number Ro is only slightly larger than 1: Ro = 1 + 6, 
where 0<(5<1. In this limit both q and p scale as 5, 
and on-site Hamiltonian ([25)1 reduces to 

H {q,p) ~ ^oqpip- q + 8) ■ (30) 

This on-site Hamiltonian describes, in WKB approxima- 
tion, a broad class of population models (that do not 
exhibit Allee effect) close to their transcritical bifurca- 
tion at 8 = [II |2(J El]. Notice that, at 5 < 1, the 
on-site dynamics exhibits critical slow-down: the char- 
acteristic on-site relaxation time becomes l/(/io<5)- As a 
result, the validity of the continuous description in space 
here demands Do 3> /j,q5: a much softer criterion than 
D > Mo- 
Let us define the characteristic diffusion length I = 
[D I '(fioS)] 1 / 2 and introduce rescaled population size Q = 
q/5, momentum P = p/S, spatial coordinate x = x/l, 
and time t = hqS t. Upon this rescaling one observes 
that the second term in the square brackets in Eq. (|24]l 
is of next order in 5 compared to the rest of terms, and 
should be neglected. The resulting Hamiltonian density 
is parameter- free, 

w = QP(P -Q + l) + Pdl Q . (31) 

Here and in the following we drop the tildes everywhere 
except in the rescaled system size L = L/l. Action (p?5)) 
becomes 

S(0) = ^s A (L), (32) 
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where 

s A (L)= dx dtP(x,t)d t Q(x,t) (33) 

JO J -oo 

is the rescaled action. The rcscalcd Hamilton's equations 
are 



d t Q 
d t P 



2QP 
2QP 



P 



Q 

p2 



d 2 x Q. 
dlP. 



(34) 
(35) 



The same WKB equations can be obtained if one ap- 
proximates, at small <5, the original master equation by 
the (functional) Fokker-Planck equation, see Appendix 
C. This is not surprising, as the validity of the Fokker- 
Planck approximation demands, in addition to > 1, 
that the probability distribution P(n, t) be a slowly vary- 
ing function of n. The latter condition boils down to 
condition \p(x, t)\ -C 1 that has been used in deriv- 
ing Eq. (|30|) . We emphasize that, far from the bifur- 
cation point, the Fokker-Planck approximation in gen- 
eral breaks down, whereas the WKB approximation still 
holds, in most of the phase space. An example is consid- 
ered in section IIVBI 

To calculate the rescaled action sa, which only de- 
pends on L, we need to find a heteroclinic connection 
AB. Deterministic steady state Q = Q s (x), P = 0, cor- 
responding to fixed point A, is given by the non-trivial 
solution of equation 



Q"(x) + Q-Q 2 = 0, 



(36) 



whereas extinction state Q = 0, P = P e (x) corresponds 
to the non-trivial solution of equation 



P"(x) +P + P 2 



0. 



(37) 



For periodic boundary conditions, see Fig. [5Jt, we ob- 
tain a;- independent solutions: Q s (x) = 1, P e (x) = — 1. 
As a result, the "extinction instanton" is x-independent 
for any system size L, and one can use the well known 
one-site WKB results [U, [H, [H]. The instanton is 
described, at any x, by the equation P — Q + 1 = 0. 
Rescaled action (l33l) along the extinction instanton is 
equal to sa = L/2. Then, using Eq. (|32|) . we find 



Hfi $T e ) ~ KS(0) 



K6 2 L _ NKS 2 
2h ~ 2 



(38) 



which is the one-site result times AT, as expected. The 
one-site result for the MTE, TjP\ is actually known with 
a higher accuracy - including a pre-exponential factor 



/<() 



ST™ 



2ir 



exp 



K5' 



Therefore, for L <C I, we obtain a more accurate result 
for T P : 



(j, ST e 



2tt 



NK5 2 



777 exp 



(39) 



as all of the system can be considered here as a single 
site. Equation holds when L < I and NKS 2 > 1. 

Now let us consider a more interesting case of absorb- 
ing boundaries: Q(0,t) = Q(L,t) = P(0,t) = P(L,t) = 
(see Appendix A). We notice that, if Eq. (f31))) has 
a nontrivial solution Qq(x), then Eq. (|37[) has a non- 
trivial solution — Qq(x). Now, Eq. (|36[) has a unique non- 
trivial solution Q s (x), corresponding to an established 
population, if L > L c = tt or, in dimensional units, 
L> L c = TrlD/inoS)} 1 / 2 , see Eq. ©. (Note that L c = nl 
here.) In this case Eq. (|37|) has a nontrivial solution 
P e (x) = —Qs(x) (fixed point B). 

Now we need to find a heteroclinic connection AB. To 
our knowledge, this cannot be done analytically for arbi- 
trary L > L c : even for relatively simple universal Hamil- 
tonian (|31[) . To solve the problem numerically, we modi- 
fied, and implemented in "Mathematica" , the algorithm 
suggested by Elgart and Kamenev [|| . The algorithm it- 
erates Eq. (|3"4")l forward in time and Eq. (|3"5)) backward in 
time. It does not involve shooting and avoids, because of 
the backward integration in time, the short-wavelength 
numerical instability caused by the presence of negative 
diffusion in Eq. (|35[) . In every iteration of Q(x,t) one 
starts, at t = 0, from Q = Q s (x) and solve Eq. ([34]) . 
with zero boundary conditions at x = and L, forward 
in time until a sufficiently long time T is reached. In 
this calculation the previous iteration for P(x, t) is used. 
Then Eq. (|35|) for P is solved backward in time starting, 
at t = T, from P = P e (x) and continuing until t = 0. 
Here the previous iteration for Q(x,t) is used, and zero 
boundary conditions at x — and L are enforced. The 
very first iteration for P is the desired final steady state 
P e (x), satisfying the zero boundary conditions at x = 
and L and corresponding to fixed point B. An example 
of numerically found instanton is shown in Fig. [SJ The 
filled circles in Fig. [10] show the numerically computed 
rescaled action s A , see Eq. (|3"3"]l . versus rescaled system 
size L/L c . 

Approximate analytic solutions are possible in two lim- 
its: L 3> L c and < L — L c <C L Cl and we will present 
now these solutions. (To remind the reader, there is no 
established population at L < L c .) 

For L ^> L c , Q s (x) and P e (x) = —Q s (x) are close 
to 1 and —1, respectively, everywhere except in bound- 
ary layers of width 0(1) at x = and L. Correspond- 
ingly, the extinction instanton is very close (up to correc- 
tions exponentially small in L) to the one-site instanton 
P — Q + 1 = everywhere except in the boundary layers. 
As a result, the rescaled action, sa — L/2 — 0(1), differs 
by a term of order unity from the corresponding result 
for periodic boundary conditions. The 0(1) correction, 
that we found numerically, is about 1.8, and its contri- 
bution to ln(/io<5T e ) is relatively large. The asymptote 
sa = ( 7T /2)(L/L C ) — 1.8 is shown in Fig. [TOJ Surprisingly, 
it works well already at quite small values of L/L c — 1. 
The next-order correction would come from the gradient 
corrections to the zero boundary conditions in space for 
q andp, see Eqs. (|A6|) . (|A8[) and (|A9|) . The expected cor- 
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FIG. 9: (color online) Numerically computed extinction in- 
stanton for scenario I (no Allee effect) close to bifurcation 
point <5 = 0. The rescaled system length is L/L c = 2. Shown, 
after 500 iterations of the Elgart-Kamenev numerical algo- 
rithm (see text), are spatial profiles of rescaled population 
size Q (a) and rescaled momentum P (b) at numerical times 
0, 3, 5, 7 and 20 (from top to bottom). The time interval 
used for iterations was < t < T with T — 50. 



criterion L ^> L c becomes stringent as the bifurcation 
point S = is approached, and L c diverges. 

Now consider the limit of < L — L c -C L c . We start 
with a perturbative calculation of Q s (x) and P e (x). Let 
e = L — L c = L — 7T<Cl. For Q s (x) we can write 

Q 8 (x) ~ ao + a\ sin x + bi cos x + 02 cos 2x , 

where a\ ~ e, whereas ao ~ 61 ~ 02 ~ e 2 . Plugging 
this ansatz into Eq. (|36[) , we obtain ao = af/2 and 02 = 
af/6. Boundary condition Q(0) = yields 61 = — 2a 2 /3. 
Now we demand Q(L) = Q(tt + e) = 0. Expanding this 
condition at small e, we obtain 



ai 



T 



L - L, 



L-TT 



(40) 



The bifurcation of the steady-state solutions, both Q s (x) 
and -P e (:z:), at L = L c = tt is, therefore, transcritical. The 
final result for Q s (x) = —P e (x), up to e 2 , is 

Q s (x)~— — h — sin x — cos x + cos 2x . (41) 
^ w 32 4 8 32 v ' 



Now let us solve perturbatively Hamilton's equations ([343 
and (|35p . Shrinking the coordinate x, 




■ x — > x . 



FIG. 10: (color online) Rescaled action, Eq. (|33[) . deter- 
mining ln(/^or e ), versus the rescaled system size L = L/L c 
for scenario I (no Allee effect) close to the bifurcation point 
5 — 0. Symbols: results obtained with the Elgart-Kamenev 
numerical algorithm (see text). Dashed line: asymptote 
s A ^ (9tt 3 /64)(L - l) 2 for < L - 1 < 1. Dotted line: 
asymptote sa — (n/2)L — 1.8 for t > 1. Inset: a blowup 
close to L = 1. 



rection to sa is O([io5 / D®) 1 / 2 <C 1. However, by virtue 
of Eq. (|32p . this correction contributes factor 0(.fr<5 2 ) to 
]n(fj,o5 T e ) . This contribution is of order of the one-site 
result for ln(/io<5 T e ) and may still be large. That is, for 
L L c the leading contribution to ln(/io<5 T e ) scales as 
L, the subleading contribution scales as L c <C i, and the 
sub-subleading contribution scales as h <C i c , the latter 
one "remembers" the lattice formulation of the problem. 
The sub-subleading correction can be calculated numeri- 
cally using the modified boundary conditions. Note that 



7T + £ 

we rewrite these equations as 



d t Q 
d t P 



2QP + Q-Q 2 



2QP - P- P 2 



2£ 
7T 
2£ 

7T 



+ 



^0,(42) 
9 2 P,(43) 



where dots denote higher order terms in e. The problem 
is now defined on the interval < x < it. We seek 
perturbative solutions in the form 

Q(x,t) = eu(x,et) + e 2 ui(x,et) + . . . , 
P(x,t) = ev(x, et) + e 2 vi(x, et) + . . . . 



In the first order in e we obtain equations 

d 2 u + u — and d 2 v + v = . 



(44) 



Their solutions, obeying zero boundary conditions at x = 
and 7r, are 



u(x, t) = <z(t) sin a; , v(x, t) = b(r) sinx . 



(45) 



where a(r) and b(r) are yet unknown functions of the 
slow time t = et. In the second order in e we obtain 



o 2 



I - — I sinx + (a 2 - 2ab) sin 2 x, (46) 

(IT 7T y 

9 2 ui+vi = - + — J sinx - (b 2 - 2ab) sin 2 x,(47) 
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subject to zero boundary conditions at x = and tt. 
The solvability conditions for Eqs. (J3H1) an d P?]) yield 
the following equations for da/dr and db/dr: 



da 2a 8 , 9 „ 

— = — - — (a 2 - 2ab) , 

db 2b 8 /l2 „ 

— = 6 2 -2a6 



(48) 
(49) 



These are Hamilton's equations for generalized coordi- 
nate a and momentum b. Hamiltonian 



H(a, b) = — ab [ b — a H — 
37T V 4 



(50) 



is of the same type as universal on-site Hamiltonian (|3U1) . 
The extinction instanton obeys b = a — 3/4, and we find 



(51) 
(52) 



4 (l + e 2 ^) ' 
3 

_ 4 (l + e- 2e */M 



One can also easily find the second-order corrections u\ 
and v\ from Eqs. (|46|) and (|47j) . but we will not present 
these formulas here. Now we calculate, in the leading 
order in e, the rescaled action ([33]) along the extinction 
instanton: 



sa 



7T+e />OC 

dx I dt p(x , t) dtq(x , t) 



e dx sin x \ drb 



da 



dr 



o 



2 J \ 4J 64 

3/4 



(53) 



This asymptote is shown in Fig. 1101 Finally, wc use 
Eqs. ([32]) and ([53]) to find the logarithm of the MTE: 



bwja^w-^g-.)-. (54) 

This result is valid, for L — L c <C L Cl when it is much 
greater than unity. This holds for sufficiently large K or 
fast migration. 



B. Extinction of particles undergoing reactions 

A 2A and 2A -> 



Here we use WKB approximation to revisit the prob- 
lem of extinction of particles A which participate in two 
on-site reactions: branching A — > 2A and annihilation 
2A 0, with rate coefficients fio and fio/K, respec- 
tively, and K ^> 1. Although still exhibiting scenario 
I of extinction, the on-site Hamiltonian for this model is 



irreducible and does not belong to the universality class 
considered in the previous subsection. In the spatially- 
independent formulation, the Fokker-Planck approxima- 
tion docs not apply for the evaluation of the MTE 
All this is because of the absence of linear decay pro- 
cess A — > (or of the linear in n small-n asymptote of 
the death rate). Here the extinction instanton, in the 
spatially-independent setting, does approach, at t — > oo, 
infinite momentum p = — oo (which, in our variables, is 
the "extinction momentum" ) . 

In spite of its degeneracy, this model is quite popu- 
lar. Its spatially-independent version was investigated in 
many papers, see e.g. Refs. [1, 0, whereas the 

spatial version was considered in Ref. @- Our objectives 
here are three-fold. First, we use this example to illus- 
trate the advantages of canonical variables Q and V (that 
arise naturally in the probability generating function for- 
malism [1, Hl[). Second, we show that, when the system 
size L is only slightly above L c , this "irreducible" model 
does reduce to the universality class considered in Sec. 
IV! Third, we use this example to compare our results 
with those of Elgart and Kamenev Q . 

We start from Eq. (|23[) for the Hamiltonian functional. 
The density w has the form of Eq. (|24|) . whereas the on- 
site Hamiltonian for processes A — > 2A and 2A — > is 
the following: 



1 



(55) 



Define characteristic diffusion length I = (D/ Hq) 1 ' 2 ■ In- 
troducing rescaled coordinate x = x/l and time t = /ioi, 
we arrive at a parameter-free Hamiltonian with density 



w = q ( e P-l) + l q 2 ( e - 2 P 
- d x q d x p + q (d x pf , 



1 



(56) 



where we have dropped the tildes. It is advantageous to 
make a canonical transformation from q and p to Q — 
qe~ p and V = e p — 1 (the shift by 1 in V preserves the 
deterministic line at "P=0) . The new Hamiltonian density 
becomes 

W = QP[V-Q + l-(l/2)QV]-d x Qd x V, (57) 

and the Hamilton's equations are 

d t Q = 2QT + Q-Q 2 -Q 2 T + d 2 Q, (58) 
d t V = OP 2 + 2QV - V - V 2 - d 2 x V . (59) 



Extinction action ([21?]) becomes 

S{0)= l -s{L), 
where L = L/l is the rescaled system size, and 



(60) 



s(L) = I dx dtV(x,t)d t Q{x,t) (61) 
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is the rescaled action. To calculate s(L) we need to find 
an instanton-like activation trajectory that exits from 
fixed point A at t = — oo and asymptotically approaches 
the proper extinction state at t = oo. Fixed point A, 
corresponding to Q = Q s (x), V = 0, is given by the 
non-trivial solution of equation 



Q"(x) + Q- Q 2 = 0, 



(62) 



with boundary conditions Q(0) = Q{L) = 0. In its turn, 
the proper extinction state Q = 0, V = V e (x) is given by 
the non-trivial solution of equation 



V"{x)+V + V 2 = 0, 



(63) 



with V(0) = V(V) = 0. Interestingly, the equations and 
boundary conditions for Q s {x) and V e {x) coincide with 
those for Q s (x) and P e {x) for the universal model of sce- 
nario I, see Eqs. (f3"6")l and (f3"T)l . In particular, equality 
V e {x) = -Q s (x) holds. 




FIG. 11: (color online) Numerically computed extinction in- 
stanton for processes A —¥ 2A and 2A — > 0. The rescaled 
system length is L/L c — 1.55. Shown, after 400 iterations 
of the Elgart-Kamenev numerical algorithm (see subsection 
II V A|) . are spatial profiles of Q (a) and V (b) at numerical 
times 0, 5, 7, 9 and 50 (from top to bottom). The time inter- 
val used for the iterations was < t < T with T = 60. 



s ~ (97r 3 /64)(L/L c - l) 2 , see Eq. (JSSJ. Surprisingly, this 
result agrees with that obtained by Elgart and Kamenev 
close to L = L c , see Eq. (45) of Ref. [f|. The reason 
for this agreement is unclear, as Elgart and Kamenev do 
not mention any boundary conditions for V at x = and 
x = L. 




FIG. 12: (color online) Rescaled action, Eq. (|61[) . versus the 
rescaled system size L — L/L c for the processes A — > 2 A 
and 2A — > 0. Symbols: results obtained with the Elgart- 
Kamenev numerical algorithm, see section llV Al Dashed line: 
asymptote s ~ (9tt 3 /64)(L - l) 2 for < L - 1 < 1. Dotted 
line: asymptote s ~ 2n(l — ln2)L — 2.36 for L> 1. 



For L^> L c functions Q s (x) and V e (x) = —Q s {x) are 
close to 1 and —1, respectively, except in the boundary 
layers at x = and L. As a result, the extinction instan- 
ton is close to the one-site instanton Q = 2(7-" — l)/(V + 2) 
everywhere except in the boundary layers, and we ob- 
tain s ~ 2tt(1 - ln2)(L/i c ) - 2.36, see Fig. [H Factor 
2(1 — In 2) comes from the solution of the one-site prob- 
lem |3, IS @i E3] ■ The leading term, proportional to L, 
coincides with that obtained by Elgart and Kamenev 0] . 
The numerically found offset 2.36 also agrees, up to 1 
percent, with their result. We can only explain this 
agreement (and the agreement in our analytic results at 
L/L c -1< 1, reported above) by assuming that Elgart 
and Kamenev did impose correct spatial boundary con- 
ditions at the edges of the system, x = and x = L. But 
then the final state V(x) in their calculations must have 
been V = V e (x) = —Q s (x), and not V = as they claim. 



POPULATION EXTINCTION: SCENARIO II, 
VERY STRONG ALLEE EFFECT 



We solved the problem numerically using the Elgart- 
Kamenev algorithm described above. An example of nu- 
merically found instanton is shown in Fig. 111! The time- 
dependent solution of this problem is different from that 
of the universal model, except when the rescaled sys- 
tem size L only slightly exceeds the rescaled critical size 
for the established population, L c = it. Here \Q\ <C 1 
and \V\ <C 1 for all < x < L and — oo < t < oo, 
and we can neglect the last term in the square brack- 
ets of Hamiltonian density (|57[) . thus arriving at univer- 
sal Hamiltonian (|31[) considered in the previous subsec- 
tion. As a result, the rescaled action in this case is again 



Here we consider the three reactions A — > and 
2A 3A. They are described, in WKB approxima- 
tion, by Hamilton's equations (f2"6")) and (|2"T|) with on-site 
Hamiltonian 

H (q,p) = Mo ( ^ + q) (e-P - 1) + ^!( e ? - 1) . (64) 

In the following we only deal with a very strong Alice 
effect in a system with periodic boundary conditions, 
see Figs. [3] and SJ Here the noise-driven population 
extinction requires a large fluctuation that creates, at 
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L > L c , critical nucleus q = q c (x). The "nuclcation in- 
stanton" (that is, a heteroclinic connection AD) can be 
found by solving Eqs. ([26]) and (f27|) with periodic bound- 
ary conditions in space for q(x,t) and p(x,t), conditions 
q(x, t — > — oo) = q-2 = 1 + S and p(x, t — > — oo) = 0, and 
conditions q(x, t — > +oo) = q c {x) and p(x, t — > +oo) = 0, 
where < S <C 1. 



A. Small and moderately large systems 

For S -C 1 the on-site deterministic dynamics is close 
to the saddle-node bifurcation. Here the unstable and 
stable fixed points, q± = 1 — d and qi = 1 + 8 are both 
close to 1, whereas the momentum p on the activation 
trajectory scales as S 2 . Expanding on-site Hamiltonian 
at small p and q — 1 we arrive at [2(| : 



H (q,p) ~ 2fi p [p + 



S 2 - Aq 2 



(65) 



where Aq = q — 1. This on-site WKB Hamiltonian, con- 
sidered already in Ref . [39| , describes a host of spatially- 
independent overdamped physical systems which exhibit 
activated escape close to a saddle-node bifurcation. Fur- 
thermore, the exact destination of the escape process 
(whether it is population extinction [20j . population ex- 
plosion [(| HH , or a switch to another metastable state 
[36j) is of no importance: it is decay of metastable state 
q = q2 which is the kinetic bottleneck of the process. 
Note that, at S -C 1, the fast-migration criterion in the 
spatial problem becomes Dq 3> hq5, as in scenario I. 

Let us define characteristic diffusion length I = 
[Df (2fioS)] 1 ^ 2 and introduce rescaled population size Q = 
Aq/5, momentum P = p/S 2 , spatial coordinate x = x/l 
and time t = 2^o5t. At 6 <C 1 the second term in the 
square brackets in Eq. ([24]) is again negligible. The result- 
ing (parameter-free) Hamiltonian density can be written 



= P[P-U'{Q) + dlQ] , 



(66) 



where U(Q) = —Q/2 + Q 3 /6 is the effective potential. 
Here and in the following we drop the tildes everywhere 
except in the rescaled system size L = L/l. Action (|29|) 
becomes 



d 3 l 
h 



(67) 



where 



s B {L)= dx dtP(x,t)d t Q(x,t) (68) 

JO J -oo 

is the rescaled action. The rescaled Hamilton's equations 
are 

d t Q(x,t) = 2P -U'{Q) + d 2 x Q, (69) 
dtP(x,t) = PU"(Q)-d 2 P. (70) 



The same WKB equations can be obtained from the 
Fokker-Planck equation, valid at small S, see Appendix 
C. 

Hamiltonian (|66)) almost coincides with the Hamilto- 
nian considered by Elgart and Kamenev [(|. The only 
difference is that the effective potential in their case was 
of the opposite sign, as they considered population ex- 
plosion rather than extinction. The procedure of finding 
the activation trajectory (heteroclinic connection AD) is 
identical in the two cases. It is based on the following 
important property of zero-energy flows in this class of 
Hamiltonians: if the quantity F(x, t) = P — U'(Q) +d 2 Q 
vanishes, at some time, for all x, then it vanishes at all 
times. This property can be easily proved by calculat- 
ing d t F(x,t) and using Eqs. (|69|) and (|70|) . Since in our 
problem Q(x,t = — oo) = 1 and P(x,t = — oo) = 0, 
equality F(x, t) — does hold. It immediately follows 
that P — dtQ-, and 



d t Q = U'(Q) - d 2 Q 



(71) 



on the activation trajectory. Equation (|71[) is a time- 
reversed version of the deterministic equation 



d t Q = -U\Q) + d 2 Q. 



(72) 



Using the relation P = dtQ in Eq. (|68|) . one obtains = 
AJF, where AF is the difference between the final (at 
t = +oo) and initial (at t = — oo) values of the rescaled 
Ginzburg-Landau free energy functional, cf. Eq. l[5]): 



F[Q(x,t)} = / dx [U(Q) + (l/2){d x Qf 



(73) 



Note that local identity F(x,t) = implies an infinite 
number of integrals of motion. Although their presence 
looks as a miracle in the WKB formalism, it is a di- 
rect consequence of integrability of the stationary Fokker- 
Planck equation in this case, see Appendix C. 

The final state, at L > L c , is the (rescaled) critical nu- 
cleus: an x-dependent solution of the steady-state equa- 
tion 



'(x) + (l/2)[l-Q 2 (x)}=0 



(74) 



subject to periodic boundary conditions with spatial pe- 
riod L. Elgart and Kamenev Q solved this equation, and 
calculated A F, in the limit of L L c . We will present 
the solution for any L > L c . The solution of Eq. (|74p . 
up to an arbitrary shift in x, can be written as 



Qc(x) = c + (b - c) sn 2 



2K(m) x 



L 



(75) 



Here b and c are two of the three real roots a(£) > b(£) > 
c(£) of the polynomial £ — £/2 + £ 3 /6 (the roots are real 
for \£\ < 1/3), sn(. . . ) is the Jacobi elliptic function, and 
K(m) is the complete elliptic integral of the first kind 
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[4C| . Furthermore, m = m{£) = (b — c)/(a 
parameter £ is determined by relation 



c), and and 



\/a — ( 



= L. 



The x-dependent solution (|75p exists at L > L c = 2tt 
or, in dimensional units, L > L c = n[2D/ ([IqS)} 1 / 2 ; this 
is what Eq. (fTTj) predicts at 5 -C 1. Note that i c = 27rZ 
here. At L > L c solution ([75)1 exhibits a single full spatial 
oscillation: its spatial period is equal to the rescaled sys- 
tem size L. At L > kL c , where k = 2, 3, . . . , this solution 
coexists with additional solutions having 2,3, ... ,k full 
oscillations. The k > 1 solutions, however, yield greater 
actions than solution ([75]) . and therefore should be ruled 
out. 

Evaluating free energy ([73")) for solution (|75|) with the 
help of "Mathematica" , we obtain sb = A J 7 = $(L/L C ) 
where function $(£) is depicted in Fig. [T3] The loga- 
rithm of the MTE is, therefore, approximately equal to 



(76) 



This result is valid when it is much greater than unity. 
This can be achieved for sufficiently large K and fast mi- 
gration. Importantly, at L > L c the MTE ceases to grow 
with system size L, so the MTE can be relatively short. 
Furthermore, for any L > L c the action spent on creating 
the critical nucleus is less than the action spent on bring- 
ing the population to the x-indepcndcnt unstable state 
q = q±. Therefore, extinction via the critical nucleus is 
(exponentially) more probable than via the state q — q±. 




FIG. 13: (color online) Function <£(£) determining the de- 
pendence of ln(/ioTe) on the rescaled system size L/L c , see 
Eq. (|76[) . for populations with a very strong Allee effect. The 
dotted line is the asymptote $(£ 3> 1) = 24/5 = 4.8. 

What are the asymptotes of this result in three charac- 
teristic regions: L < L c , < L — L c -C L c and L S> L c l 
Instead of using asymptotics of the elliptic functions, one 
can directly solve Eq. (|74"|) in these regions. At L < L c 
the critical nucleus gives way to the ^-independent un- 
stable solution Q = —1. Here we obtain AJ 7 = 2L/3, 



2LK5 3 2NKS 3 



3/i 



(77) 



which is the one-site result |6|, [2Cj times TV as expected. 
Again, for L -C L c we can use a more accurate one-site 
result [20| and obtain 



fi Q 8T e ~ 7T exp ( ^NKS 3 



(78) 



as the whole system can be considered as a single site. 
This result is valid when L«L C and NK6 3 > 1. 

At L = L c a weakly inhomogeneous critical nucleus 
emerges via a super-critical bifuraction. At < L — L c <C 
L c the critical nucleus, in the rescaled variables, is 



,'2nx\ A 2 A 2 

-fx) ~ —1 + A cos —=— H — cos , 

L J 4 12 \ L J 



( 47Tx\ 



where 



(79) 
(80) 



Here we obtain 



18 / L 

'-ii 



and so 



hi(/x 6 T e 



2NK5 3 



18 



L 



(81) 



Finally, at L 3> L c the critical nucleus can be approxi- 
mated by its asymptote at L — > oo: 



3 (x) ~ 1 - 3coslT 2 (x/2) 



(82) 



In this limit, mathematically identical to the one consid- 
ered by Elgart and Kamenev @ in the context of popu- 
lation explosion, we obtain sb = AJ 7 = 24/5, and so 



ln(/x 5 T e ) 



24KS 3 l \2K8 3 L r 



■5// 



hirh 



(83) 



Note that criterion L ^> L c becomes stringent as the 
bifurcation point S = is approached, and L c diverges. 

Now we see that the asymptotes of $(£) are the fol- 
lowing: 



3 S : 




, 0<f-l«l, (84) 

e»i. 



Before concluding this subsection we note that, al- 
though we succeeded in calculating action sb analyti- 
cally, the calculation of the nucleation instanton, Q(x,t) 
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FIG. 14: (color online) A numerically computed nucleation in- 
stanton (heteroclinic connection AD) for scenario II (a strong 
Allee effect) close to bifurcation point 5 = 0. The rescaled 
system length is L/L c ~ 1.54. Shown are spatial profiles of 
the rescaled population size Q at numerical times 0, 11, 12, 
13, 14 and 20 (from top to bottom). The last profile is the 
critical nucleus for this system size. 



and P(x,t), demands (rather simple) numerics. Q(x,t) 
is described by the time-reversed version of reaction- 
diffusion equation (|72[) . Therefore, for a given system 
size, one can solve Eq. (|72[) with periodic boundary condi- 
tions numerically, starting from the critical nucleus (with 
a slight positive offset) and advancing the solution until 
it converges close to Q = 1. The instanton solution is 
then readily obtained via time reversal, whereas P(x, t) 
can be found from P(x,t) = dtQ{x,t). An example of 
instanton, found in this way, is depicted in Fig. [T4l 



B. Very large systems: single versus multiple 
nucleation 

In the previous subsection we evaluated ln(/io<5T e ) un- 
der condition that the probability of creating more than 
one critical nucleus during the traverse time of the deter- 
ministic extinction fronts through the population is neg- 
ligible. The rest of parameters being fixed, this condition 
is always satisfied at sufficiently large K. If we instead 
fix K 3> 1 and increase L, we will arrive at the regime 
when additional critical nuclei typically appear while the 
extinction fronts still run through the population. The 
present subsection deals with this regime. Importantly, 
the assumption of quasi-stationarity, sec Eqs. (|15[) and 
(|16p . does not hold in this regime, except for the purpose 
of calculation of the rate of formation of a single critical 
nucleus in the phase q% = 1 + 5. The latter is given, at 
L ^> L c and t >• l/(p S), by Eq. ([55]) that we rewrite 
here, up to pre-exponential factors, as 



1 

t1 



cc p S e KSo ' where So 



12S A L C 



(85) 



As l/T e is exponentially small, the nucleation acts are 
rare, and we can assume that they are statistically in- 
dependent and Poisson-distributcd. This implies that, 
in sufficiently large systems, L 3> £ c , nucleation rate 



(1851) includes a pre-exponential factor proportional to L 
(that we did not care about previously but must account 
for now). Correspondingly, the nucleation rate density p 
(that is, the nucleation rate per unit length of the system) 
is independent of i, and we can represent it as 



p 5 



where 



= Re 



-KS 



< 1 



(86) 



and R is a dimensionless prc-factor that depends on di- 
mensionless parameters K, 5 and ho/Dq. Importantly, 
the nucleation problem is mathematically equivalent to 
the over-damped limit of theory of homogeneous nucle- 
ation due to Langer [27}, see Appendix C. The theory 
of Langer corroborates our argument that, for L L c , 
the nucleation rate is proportional to L. Furthermore, 
his theory makes possible to calculate pre-factor R ex- 
plicitly. We will not need the pre-factor, however, as we 
are only interested in the leading-order approximation 
for ln(p, 5T e ). 

With nucleation rate density (|86p at hand, we now con- 
sider the following problem. Let at t = the whole sys- 
tem, with L 3> L c , be in the populated state qi = 1 + S. 
The nucleation rate density p is independent of x and 
t. After a critical nucleus (of size ~ L c <C L) develops, 
two deterministic extinction fronts form and propagate 
in both directions with speed c ~ y/ pqD/2, see Eq. (|12l) . 
What is the probability V q2 (xq, to) to still observe q = qi 
at point x = xo at time t = to > 0? For this to hap- 
pen, no critical nucleus should have appeared in space- 
time domain G within the event horizon produced by 
the two incoming extinction fronts. Taking into account 
the finite size L of the system and the periodic bound- 
ary conditions, wc find that, for 2cto < L, domain G is 
determined by conditions 

\x - x \ < c(t - t) and < t < to , (87) 

whereas for 2cto > L it is determined by conditions 

< x < L for < t < t - — , and 



2c 



L 



\x — x \ < c(t — t) for t - — <t<ta, 

2c 

see Fig. [15] Indeed, because of the periodic boundary 
conditions (which bring translational invariance of the 
problem) we can always choose xo = L/2, so V q . 2 (xo,to) 
is actually independent of xo- The space-time area a (to) 
of domain G is equal to 



ctr, 



for 2c£o < L . 



By virtue of the Poisson statistics, we obtain 



v (t)-e-^-r xp y pct) T, 

/?2W " | exp(-pLt+^r) for 2ct > L 



for 2ct < L , 
L. 
(90) 
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FIG. 15: Event horizon produced by two extinction fronts: at 
2ct < L (a) and 2ct > L (b). 

The MTE can be calculated from T e = V q2 (t) dt, 
and we obtain 

r -i/5-'6/s <»> 

where erf(. . . ) is the error function [ioj . Using expres- 
sion (f8"rJ)) for p, we can rewrite Eq. (pTTj) as 

where the characteristic length scale L* (which is expo- 
nentially large in K) is defined as 

L, = ./X j-i/4 Lc = 25- 3 /4 . (93) 

V 7T/0* V MOP* 

For i c « L « L, the first term in Eq. (|92[) can be ne- 
glected, and we recover, up to a pre-exponent, the result 
from subsection fVl A: poST e ~ [L c /L)e °. However, for 
exponentially large systems, L ^>> the second term 
is negligible, whereas erf(i/L*) — > 1, and we arrive at 
Po8T e ~ l/y/fa ~ e KS °/ 2 . Note that, if we interpret this 
new asymptote as exponential of some effective WKB 
action, this action will be twice as small as the action 
obtained for L< L,. A sketch of the overall dependence 
of the MTE on the system length is presented, on a log- 
log scale, in Fig. [Tf)J Evident is a maximum of the MTE 
at L much larger than L c but much smaller than L*. 

The exponentially large characteristic length scale L* 
comes, quite naturally, from the balance between the nu- 
cleation rate, pL = po5p*L/ 'L c and the traverse rate of 
the deterministic extinction fronts through the system, 
~ c/L. At L < L, the population is typically going ex- 
tinct via formation of only one critical nucleus, whereas 
at L L„ multiple nucleation acts typically occur. 



Finally, the behavior of the MTE versus L, depicted in 
Fig. [T51 can be understood as follows. At L c < L < L» 
it is the formation of a single nucleus that serves as a 
bottleneck of the extinction process. As the formation 
rate of the nucleus goes up linearly with L, the logarithm 
of the MTE goes down linearly with L in this regime. 
The linear decrease reaches a plateau at L > L* when 
multiple nucleation acts occur, and multiple extinction 
fronts are at work. 

th 6T e I 1 




L = h L = L C L = L 



FIG. 16: Shown, on a log-log scale, is a sketch of system- 
size dependence (|92|l of the rescaled MTE of a population 
exhibiting a very strong Allee effect. 



VI. DISCUSSION 

When an isolated stochastic population resides in a 
refuge of a large but finite size, it ultimately goes ex- 
tinct with certainty. We have developed WKB approx- 
imation to the quasi-stationary multi-variate probabil- 
ity distribution of the population sizes and arrived at an 
effective Hamiltonian mechanics that encodes the most 
probable path the population takes on the way to ex- 
tinction, and enables one to evaluate the mean time to 
extinction (MTE). The most general, spatially discrete 
version of WKB equations employs lattice Hamiltonian 
([2T]) and is valid for (almost) any relation between the 
migration rate coefficient Dq and the characteristic rate 
coefficient po of the on-site birth-death dynamics. For 
example, one can use these equations to address an inter- 
esting regime, in populations with an Alice effect, where 
discreteness of the lattice and a low migration rate con- 
spire to cause propagation failure of deterministic fronts 
(of either extinction, or colonization), see Ref. [4l| and 
references therein. If the migration is much faster than 
the on-site population dynamics, it can be described as 
diffusion, and one arrives at an effective continuous clas- 
sical mechanics, Eqs. (|23[) - (p3|) . where one has to find 
an activation trajectory: the most probable path of the 
population to extinction. In the absence of Allee effect 
(extinction scenario I) and for a very strong Allee effect 
the most probable path to extinction is an instanton - 
a proper heteroclinic connection in the functional phase 
space of the system. 

The extinction dynamics, and the MTE, can be very 
different depending on whether the population exhibits, 
or not, Alice effect, as well as on the conditions at the 
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refuge boundaries. The most dramatic differences appear 
for a sufficiently large system size, L> L c . In this case, 
in the absence of Allee effect (extinction scenario I), the 
MTE continues to grow exponentially with the system 
size. When a very strong Allee effect is present, however, 
extinction occurs via formation of a critical nucleus, and 
the MTE becomes, up to a pre-exponcnt, independent 
of the system size. We have obtained detailed results by 
assuming that the birth and death rate coefficients are 
such that the system is close to its characteristic bifur- 
cation (transcritical or saddle-node in scenarios I and II, 
respectively). In these cases one obtains universal Hamil- 
tonians (joTj) and (|66|) . describing two broad classes of 
population models: without Allee effect, and with a very 
strong Allee effect, respectively. We have also revisited 
the model system A —> 2 A and 2 A and shown that, 
close to the critical system size L c , this system belongs 
to the universality class described by Hamiltonian pip . 

For a very strong Allee effect we have mapped the 
extinction problem into the over-damped limit of the- 
ory of homogeneous nucleation due to Langer, where the 
corresponding stationary Fokkcr-Planck equation is inte- 
grable. This connection gives a natural explanation to 
the intcgrability of the zero-energy WKB problem con- 
sidered in section [V] In very large systems the MTE 
starts to go down with the system size, so there is an 
optimal refuge size for which the MTE is maximum. At 
still larger systems the dependence of the MTE on the 
system size reaches a plateau. Here multiple nucleation 
acts occur, and multiple extinction fronts arc at work. 

For extinction scenario I, the Elgart-Kamenev algo- 
rithm of forward and backward iterations @ yields ac- 
curate results for the MTE, and for the most probable 
path to extinction. An efficient algorithm that would deal 
with spatial populations exhibiting extinction scenario II 
is unavailable as of present. This hinders progress of 
theory beyond the completely integrable case of a very 
strong Alice effect. The weak- Alice-effect regime remains 
terra incognita. This includes evaluation of the MTE for 
periodic or reflecting boundaries, see scctionlH"lC2, where 
the choice between different possible paths to extinction 
is not obvious. 

In the general part of our derivation, section IIII1 we 
presented the WKB theory for single-step birth-death 
processes, and in one spatial dimension. A generaliza- 
tion to multiple-step processes (such as a simultaneous 
birth or death of more than one individual) is straight- 
forward, see section IIV B[ and was already introduced 
in Ref. [Q. Higher spatial dimensions can be also taken 
care of. For scenario I this was observed in Ref. [6f . For 



scenario II (a very strong Alice effect), the problem re- 
mains integrable in higher dimensions, except that the 
critical nucleus must in general be found numerically. 
More challenging generalizations include multiple popu- 
lations (competition, predation, infection/recovery, etc.) 
and environmental noise. 

Our WKB calculations were based on the assump- 
tion that the classical action for the one-site problem is 
much greater than unity. This assumption necessitates 
K 3> 1. Our fast- migration results, however, strongly 
suggest that this criterion can be relaxed. For example, 
it is obvious that, for homogeneous-in-space regimes of 
extinction, one can treat the whole system as a single 
site, see Eqs. ([59")) and ([75)1. and it is the resulting ac- 
tion for the whole system that only needs to be large for 
WKB theory to hold. As N ^> 1, the latter condition 
can be satisfied even for K < 1. For inhomogeneous- 
in-space extinction regimes, it should suffice to demand 
that the action contributed by regions whose spatial di- 
mension is of order of the characteristic diffusion length 
I ~ L c be much greater then unity. This necessitates 
KL c /h ~ K(D / /io) 1 / 2 > 1. For a fast migration, 
Do ^> fio, this condition is much softer than K ^> 1, 
and it may be be further relaxed close to characteristic 
bifurcations of the on-site Hamiltonians. 

Put in a more general context, this work dealt with 
rare large fluctuations in spatial stochastic systems far 
from thermal equilibrium. The last decade has seen a 
surge of interest in a similar class of problems in the con- 
text of steady-state currents in spatial systems of inter- 
acting particles, driven by reservoirs at the boundaries, 
see e.g. Ref. [42( and references therein. WKB approxi- 
mation, bringing about the Hamilton- Jacobi or, alterna- 
tively, Hamilton's formalism in a functional phase space, 
has been instrumental in the analysis of those systems as 
well HHH. 
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Appendix A. Absorbing boundaries: governing 
equations and spatial boundary conditions 

Here we consider a refuge with absorbing boundaries. 
The individuals can exit the refuge through its edges i = 
1 and i = N (with the same migration rate coefficient 
Dq as in the bulk), but no individuals can enter the sites 
i = 1 and i = N from outside. In particular, this setting 
models the extreme situation when the conditions outside 
of the refuge are so harsh that the individuals die there 
instantaneously. In this case master equation (fT3]l needs 
to be replaced by the following one: 
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d t P(n, t) = Yj X ( n * - 1 ) P ( A ' n * -!»*) + + l)P(n, + 1, t) - [A(n») + fx(m)]P(n, t) 

i=l 

W-l 

+ L»o 51 K-i + l)P(n, Hi-i + l,m - 1, i) + (m+i + l)P(n, - 1, n i+1 + 1, t) - 2niP(n, i) 
+ L>o + l)P(n, ni + 1, t) + (n 2 + l)P(n, ni - 1, n 2 + 1, i) - 2mP(n, i) 

+ (riiv-i + l)P(n, njv-i + 1, % - 1, i) + (njy + l)P(n, % + 1, t) - 2riArP(n, f)] . (Al) 



Going over to the eigenvalue problem, as in Eq. (fTT|). 
and applying WKB approximation (fT9")l , we obtain the 



This lattice Hamiltonian \cf. Eq. (|2ip] holds for any re- 
lation between Do and /^o [331 ]. Now let us consider the 
limit of Dq [j,q. Here for smooth solutions such as, e.g. 
activation trajectories, one has \pi — <C 1. Proceed- 
ing as in Sec. IIII1 we can Taylor-expand the migration 
part of the Hamiltonian: 

jv 

H m {q, p) = D Y [ ~ (li - Qi-i) (Pi - Pi-i) 

+D [qi (e- pi -l)+q N ((T>« - l)] .(A3) 
Now consider the Hamilton's equation for dpi/dt: 

^± = -D (e-^-l+p 2 - Pl )+... , (A4) 

where . . . denote small corrections coming from the on- 
site Hamiltonian O(fio) and higher-order terms in p 2 — 
Pi . The characteristic time scale of the dynamics of the 
system (for example, on an activation trajectory) is /Uq 1 
(or longer when a bifuraction is approached). Therefore, 
the left hand side is small, and we obtain, in the leading 
order in fio/Do, 

e- pi -i +P2 _ Pl ~o. (A5) 

As \ P 2 — pi I -C 1, the only way to satisfy this con- 
dition is to assume that p\ <C 1 which yields, up to 
small corrections, p 2 = 2pi. Rewriting this relation as 
Pi ~ (j>2 ~ Pi) =0 and going over to continuous descrip- 
tion, we obtain 

p{x = 0,t)-h d x p(x = 0, t) = (A6) 



I 

following WKB-Hamiltonian: 



-1). 

(A2) 

I 

or, in the leading order, simply p(x = 0, t) = 0. 

Now wc consider the Hamilton's equation for dqi/dt. 
Up to small corrections, wc obtain 

^ = Do (q 2 - 2qi) + ..., (A7) 

so again q 2 = 2qi+ small corrections. This yields, in the 
continuous description, 

q(x = 0,t)-hd x q(x = 0,t) = (A8) 

or, in the leading order, q(x = 0, t) = 0. Repeating these 
arguments for site i = N we obtain 

q{L,t) + hd x q(L 7 t) = 0, p(L,t)+hd x p(L,t) = (A9) 

or, in the leading order, q(x = L,t) = p(x = L,t) = 0. 

When going over to continuous description in the bulk, 
one arrives at the same continuous Hamiltonian (| 2 3 1) - ([25]) 
as in the periodic case. 

We note that the gradient terms that appear in 
Eqs. (|A6|) . (|A8[) and (|A9|) can be legitimately taken into 
account as small corrections to the zero boundary con- 
ditions for q and p. Indeed, they are of relative order 
(Ho/Do) 1 / 2 [because the characteristic length scale of the 
problem is I ~ (D/fio) 1 ^ 2 = h(Do/ Ho) 1 ^ 2 ]: whereas the 
omitted terms - both in the boundary conditions and in 
the Hamilton's equations in the bulk - are much smaller: 
of order /xo/Do- Close to the transcritical bifurcation, see 
Sec. IIV1 one should replace /zo by fio$ hi these estimates. 

Zero boundary conditions for the momentum also ap- 
pear in the context of large deviations in open systems of 
interacting particles, driven by reservoirs at the bound- 
aries [Mllg. 



N N-l 

P"(q, P) = Mo E [%0 ( ePl " !) + ( e ~ P< - !)] + D ° E [ 



Qi-i e" 



-1 



-Pi+i 



+ Do [qi {er^ -l)+q 2 (e^ 2 - l)] + D [q N _i (f»-P»-* - l) + q N ( e "P" - l)] 
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Appendix B 

The statement that the activation trajectory in extinc- 
tion scenario I must be a heteroclinic connection AB in 
the functional phase space {q(x),p(x)} relies on the pres- 
ence and linear stability properties of fixed points - that 
is, steady-state solutions with specified boundary condi- 
tions in space of Eqs. (|26|) and (|2"T|) . These are described 
in Appendix Bl. This information is then used in Ap- 
pendix B2 to prove the statement. Appendix B3 presents 
a linear stability analysis of fixed points A, C and D of 
scenario II. 



1. Scenario I: Functional fixed points and their 
linear stability 

As a typical example of scenario I, we consider univer- 
sal Hamiltonian (|31[l introduced in section IIV Al There 
are three zero-energy fixed points here. 



(i) Fixed point A: Q 



,(x), P = 



Here we put Q(x,t) = Q s (x) + q(x,t) and P(x, t) = 
p(x, i) and linearize rescaled Eqs. (fM)) and ([33|) with 
respect to q and p. The linearized equations are 

d t q = -2Q s (x)q + q + dlq + 2Q s {x)p, (Bl) 
d t p = 2Q s {x)p-p-dlp, (B2) 

subject to zero boundary conditions for q and p at x — 
and x = L > L c . Start with Eq. (|B2[) and look for eigen- 

Eigenfunctions 

, (B3) 



modes of the form p(x,t) = e ip{ x ) 
ip(x) satisfy Schrodingcr equation 

i)"(x) + [<£- V(x)]ip(x) -- 



where <£ = E + 1, and V(x) = 2Q s (x). 

Now we will prove a simple comparison theorem. Con- 
sider an auxiliary equation: 



V>"0) + 



<B-\v{x) 



V>(x) = . 



(B4) 



By virtue of Eq. (|36j). it has nontrivial solution ip(x) oc 
Q s (x) at (E = 1. This solution obeys zero boundary con- 
ditions at x = and L and has no nodes inside the in- 
terval < x < L. As a result, £ = 1 is the lowest 
eigenvalue of the auxiliary problem. Now, our original 
potential V(x) in Eq. (|B3I) is higher everywhere, except 
at points x = and L, than auxiliary potential V{x)/2. 
Therefore, the lowest eigenvalue of original problem (|B3[) 
is strictly greater than the lowest eigenvalue of the aux- 
iliary problem. Hence min (£ > 1, and so all eigenvalues 
E are positive. 

Now we turn to Eq. (|B1[) . Here q(x,t) is forced by 
term 2Q s (c) p(x, t). Let us expand both the forcing and 



the solution that we are seeking in the complete set of 
eigenfunctions i/j n (x) of the momentum: 

oo 

2Q S {x)p(x, t) = J2 bne E -^n (x) , (B5) 



and 



9(M) = ^2 fn(t)tp n (x) , 



(B6) 



n=l 



We obtain equation 



-TT + E nfn = 

at 

for fn(t), whose general solution is 



E n t 



f n {t) = me"*"' + 



E n t , „E n t 



2E n 



As E n > for all n = 1, 2, . . . , we see that deterministic 
hyper-plane p = is the stable manifold of fixed point A, 
as expected from the deterministic theory. The unstable 
manifold involves a non-zero p. 

(ii) Fixed point B: Q = 0, P — P e (x) 

Here we put Q = q{x,t) and P = P e {x) +p(x,t). The 
linearized equations are 

8 t q = -2Q,(x)q + q + %q, (B7) 
d t p = 2Q s {x)p-p-d 2 x p-2Q s {x)q. (B8) 



where we have used the fact that, for Hamiltonian (|3ip . 
P e (x) = —Q s (x). The analysis here is very similar to 
that for fixed point A. We first consider Eq. (|B7|) and 
look for eigenmodes q(x,t) = e Et <fi(x). We observe that 
Eq. (|E>7|) for q coincides, up to a sign, with Eq. (|B2[) for 
p. As a result, eigenvalues E n are "mirror images" of E n , 
considered in the context of fixed point A and, therefore, 
all of them are negative. 

The analysis of forced Eq. (|B8[) for p closely follows 
that for forced Eq. (|Blj) for q. Wc expand the forcing 
and the solution in the complete set of eigenfunctions 
(j> n (x) of q: 



2Q s (x)q(x, t) = - d n t 



-E n t 



4>n{x) 



and 



P(x,t) = 22 9n{t)<t>n{x) , 



(B9) 



(BIO) 



and obtain 



dg n 
dt 



-Ent 



23 



The general solution is 

As E„ > for all n = 1,2,..., hyper-plane q = is 
the unstable manifold of fixed point B, whereas its stable 
manifold involves q(x) > 0. 

(Hi) Fixed point C: Q = P = 

Here the linearized equations are 

d t q = q + 8 2 x q, (Bll) 
d t p = -p-Olp, (B12) 

and the eigcnmodes are elementary: 

q = Ae r "i*sin— , (B13) 
L 

p = Be 7 "2* sin^=p . (B14) 

where T ni = 1 — nfn 2 /L 2 , j n . 2 = — 1 + n 2 Tr 2 /L 2 , and 
n\,n2 = 1,2,.... We are only interested in solutions 
with q{x, t) > 0, therefore the only allowed mode for q is 
the fundamental: m = 1. This mode, at L > i c = 7r, 
is unstable. As a result, region q(x) > in a vicinity 
of fixed point C belongs to the unstable manifold of this 
fixed point. 



Fixed point B has unstable manifold q(x) = 0, and 
a stable manifold that belongs to domain q(x) > 0. 
Each of the two manifolds is A^-dimensional in the lat- 
tice formulation and is embedded into zero-energy hyper- 
surface H{q(x),p(x)} = 0. Now we see that we need to 
find a trajectory going from fixed point A to fixed point 
B. This trajectory must belong to both hyper-surfaces 
£a and E#. In the lattice formulation, each of the two 
hyper-surfaces is iV-dimcnsional [and is embedded into 
(2N — l)-dimensional hyper-surface H{q(x),p(x)} = 0]. 
Therefore, hyper-surfaces T,a and can intersect, in 
general, only along a finite set of one- dimensional curves 
which are trajectories generated by Eqs. ([26]) and (j27|) . 
These are heteroclinic connections. If there are more 
than one such connections, the one with the minimum 
action along it determines the MTE. 

We observed numerically that there is exactly one het- 
eroclinic connection AB in two different examples: for 
universal Hamiltonian pip and for the set of reactions 
A — > 2 A and 2 A — > 0, see section HVl This property ap- 
parently holds in a broad class of single-population sys- 
tems exhibiting extinction scenario I. 

3. Scenario II: fixed points A, C and D and their 
linear stability 

In the limit of a very strong Allee effect we only need to 
investigate the linear stability properties of fixed points 
A, C and and D. 



2. Scenario I: Activation trajectories are 
heteroclinic connections 

Consider extinction scenario I in a system with absorb- 
ing boundaries at x = and x = L > L c . As one can 
see from Eqs. (|2"rJ|) and ([2"T]) . hyper-planes q(x) = and 
p(x) = are invariant manifolds. Each of them is embed- 
ded into zero-energy hyper-surface H{q(x),p(x)} = 0. 
Therefore, hyper-plane q(x) = cannot be reached from 
domain q(x) > except via fixed points of Eqs. (|26[) 
and (|2T|) that belong to hyper-plane q{x) = (and have 
finite p), or alternatively via p(x,t) = — oo. There are 
exactly two fixed points belonging to hyper-plane q = 0: 
B and C. As domain q > in a small vicinity of fixed 
point C belongs to its unstable manifold, fixed point C is 
unreachable from domain q(x) > 0. On the contrary, the 
stable manifold of fixed point B does include q(x) > 0. 
Therefore, a trajectory can exist that asymptotically ap- 
proaches fixed point B at t — > +oo. 

Now consider trajectories of Eqs. (|2l)|) and (f2~T)) that 
come into hyper-plane q(x) = [at a finite time, and 
simultaneously at all points of the open interval (0, L)} at 
p{x,t) = — oo. One can show that H{q(x 7 t),p(x,t)} > 
for such trajectories, and so they cannot start from fixed 
point A. This is similar to what happens in spatially- 
independent but multi-population systems (22j . 



(i) Fixed points A and C: q = qi or 0, p = 0. 

Linearizing Eqs. (|2f)|) and (j2~T)) around fixed point A, 
we obtain 

d t Sq = ^f'(q 2 )8q + Dd 2 x 8 qi (B15) 
dtp = -lx f(q 2 )p-Dd 2 x p. (B16) 

As /'((ft) < 0, one can see that deterministic hyper-plane 
p = is a stable manifold of fixed point A, whereas p ^ 
is an unstable manifold. The same results hold for fixed 
point C. 

(ii) Fixed point D: q — q c (x + const), p — 0. 

For periodic boundaries, there is a one-parameter fam- 
ily of fixed points D corresponding to an arbitrary shift 
with respect to x. Because of this degeneracy, there are 
two eigenmodes that correspond to zero eigenvalue. One 
of them is neutrally stable: 

q(x) - q c (x) = const q' c (x) , (B17) 
P = 0; (B18) 

it corresponds to an infinitesimal shift in x of the crit- 
ical nucleus q = q c (x). The other one is algebraically 
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unstable, as it grows linearly in time: 

q(x)-q c {x) = aq' c (x)t + ip(x), (B19) 
P = Cq' c (x), (B20) 

where ip{ x ) obeys the periodic boundary conditions. We 
skip here the exact form of function ip(x), as well as the 
expression for non-zero constant a. 

As a result, in the lattice formulation we would have 
an (N — l)-dimensional stable manifold T,u that con- 
tains point D and leaves hyper-plane p{x) = 0; a 
one-dimensional neutrally stable manifold, belonging to 
hyper-plane p(x) — 0, and an TV-dimensional unstable 
manifold with one of its tangent vectors belonging to hy- 
perplane p(x) = 0. The neutral manifold corresponds to 
a one-dimensional line of fixed points D, parameterized 
by the exact position of critical nucleus on interval (0, L). 



tion [32j . Unfortunately, it can only give accurate results 
for the MTE when the system is sufficiently close to bifur- 
cations describing emergence of established populations 
(for spatially-independent problems this was observed in 
Ref. [7[|). Indeed, only in this case the probability distri- 
bution of the population size is a slow varying function 
of the populations size at all relevant population sizes, so 
that the (truncated) van Kampen system-size expansion 
[32^ becomes accurate. Here we derive (lattice versions 
of) Fokker-Planck equations close to the characteristic bi- 
furcations of extinction scenario I and, for a strong Allee 
effect, scenario II. These systems are analyzed, in WKB 
approximation, in sections IIV Al and |V] respectively. We 
also point out to mathematical equivalence between the 
problem of population extinction for a very strong Alice 
effect and the over-damped limit of theory of homoge- 
neous nucleation due to Langer 27 (. 



Appendix C. Fokker-Planck equations for scenarios I 

and II 

The Fokker-Planck approximation is a commonly used 
large-population-size approximation to the master equa- 



We start from a formal truncated Taylor expansion 
of the multivariate probability distribution in time- 
dependent master equation (|13[) . The migration terms 
in Eq. (|T3|) become 
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Here we have assumed, for concreteness, periodic bound- crything together. The result is a formal Fokker-Planck 
ary conditions. Now we also expand the on-site terms, equation 
employ Eq. ([2]) , go over from rij to q± = rii / K and put ev- 
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Even for K 1, this equation is only valid, in general, 
around the attracting fixed point of the deterministic rate 
equations, where the long-lived population distribution 
resides. In addition, this equation can hold, for extinc- 
tion scenario II, in the region around the repelling fixed 



J 



point. Importantly, it does become accurate for all popu- 
lation sizes rii ^> 1 when the system is sufficiently close to 
a bifurcation corresponding to emergence of established 
populations. For the SIS model (scenario I) we can as- 
sume <5 <gC 1 and obtain, after some algebra, 



dP(q,t) 
dt 



Mo 
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where we have rewritten the first-derivative migration Similarly, for three reactions A — > and 2A <=± 3A (sce- 
term in a divergence form and neglected the second- nario II), we obtain, for 5 <C 1 (a very strong Alice effect): 
derivative migration term (as it is of next order in S). 



dt 



N „ 

^ E {dq- [(* " 1 + " 1 " 6)P(q, t)] 



2 <3 2 P(q,t) 
K dq? . 



N 



~ D ° E a~~ ~ 2 * + 



(C4) 



r 



The neglected second-derivative migration term includes 
an additional S 2 factor. In the fast-migration limit, one 
can replace the lattice formulation by a continuous one, 
arriving at functional Fokkcr-Planck equations close to 
the bifurcations of scenarios I and II. 

Each of Eqs. (|C3[) and (|C3[) can be rewritten as a 
continuity equation, 



dP 
~dt 
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For scenario I (at 5 <C 1) the probability flux is 

^-nr^r-Tir^' (C6) 

c% K dqi 

with free energy 
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For scenario II we have, also at 5 <C 1, 

aF(q) p 2/io 9 p 



(C7) 
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Fokker-Planck Eq. (|C8|) is closely related to the Fokker- 
Planck equation that appears in the homogeneous nucle- 
ation theory of Langer, see Eqs. (2.15) and (2.16) of his 
paper (27j . This close relation turns into a full equiv- 
alence if one goes, in Langer's equations, to the over- 
damped limit, Aij = 0, sets T = kT = 2/io/ A, and 
specifies free energy .F(q) as in our Eq. (|C9|) . It is cru- 
cial that the multi-dimensional effective force, that ap- 
pears in probability flux (|C8|) . is potential, whereas the 
diffusion coefficient 2[1q/K is q- independent, as if com- 
ing from additive white Gaussian noise in the equivalent 
Langevin formulation of the problem. In this case the 
multi-dimensional Fokker-Planck equation is integrable 
for the purpose of calculating the stationary distribu- 
tion, and this integrability is closely related to detailed 
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balance property, sec e.g. Rcf. 
gration one can go to the continuous limit and rewrite 
free energy (|C9|) . upon rescaling, as our Eq. (|T3"|) . The 
corresponding infinite-dimensional problem remains inte- 
grable for the purpose of finding the stationary solution of 
the (functional) Fokker-Planck equation. More precisely, 
the conservation law F{x,t) = 0, that appears in section 
IV Al immediately follows, in WKB approximation, from 
the continuous version of zero probability flux condition 
Ji = that solves Eq. (|C5|) . This clarifies the reason 
behind integrability of the zero-energy WKB problem, 
considered in section IV A i m the context of a very strong 
Allee effect (and, in Ref. [g, in the context of population 
explosion). 

The situation is less fortunate for extinction scenario 
I. Although the multi-dimensional force, entering prob- 
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ability flux (|C6|) . remains potential, the diffusion coeffi- 
cient here is g-dependent, as if coming from multiplica- 
tive white Gaussian noise in the Langcvin formulation. 
As a result, the problem of finding the stationary distri- 
bution is non-integrable here, and this manifests itself in 
the non-integrability of the zero-energy WKB problem 
considered in section HV Al 



The formal equivalence between Eqs. (|C5[) . (|C8|) and 
(|C9|) for a very strong Allee effect and the equations of 
theory of homogeneous nucleation makes it possible to 
go beyond the leading WKB-order and calculate the sub- 
leading correction (a pre-exponential factor) to the nu- 
cleation rate 271. 



